Homing missile guidance and estimation algorithms against advanced maneuvering targets

ABSTRACT

A system, processing structure and method for performing vehicle intercept includes target acceleration estimation or blocking via a filter prior to deriving an interceptor acceleration time history via a guidance law. In the processing structure where the unknown target accelerations are blocked from affecting estimates, a guidance law is generated, stored and retrieved accounting for worst case target acceleration. In the processing structure where the target acceleration is presumed constant with a constant bank angle, the guidance law is closed form, exploiting constant heading angular velocity, and thus generated online.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of provisional application No. 60/696,739, filed Jul. 5, 2005, the disclosure and appendices of which are hereby incorporated by reference herein, in their entirety, for all purposes.

BACKGROUND OF THE INVENTION

1. Technical Field

The present invention, in its several embodiments, relates to vehicle intercept methods and intercept processing structures and more particularly to the methods of, and processing systems for intercept estimation and interceptor guidance applied to maneuvering targets.

2. State of the Art

Because the process of proportional navigation of an interceptor does not yield a solution to the intercept of maneuvering targets, several advanced guidance laws and intercept geometry filtering techniques have been developed in attempts to accommodate the targets having time variation in the magnitude or direction of their respective velocity vectors. One approach in the construction of the filter and guidance law that has been applied to maneuvering target has been to use a target model, typically a discrete mathematical representation, in which the target is presumed to perform particular maneuvers. Based on such a target model, a guidance filter can be developed to estimate the relative target-to-interceptor position and velocity and target acceleration, as state estimates for example. Furthermore, the associated guidance law, which generates the interceptor's acceleration command using the state estimates, can be developed subject to the presumed engagement dynamics. However, a considerable drawback wit this target model approach is that the target acceleration cannot be easily modeled because the time history of the target acceleration is inherently a jump process whereby the acceleration levels and switching times are unknown a priori. Therefore, there remains a need for a more realistic target model and the associated filter and guidance law based on this target model.

Instead of presuming a target model, another approach for constructing the guidance law is to presume the worst possible target acceleration within the interceptor engagement scenario; i.e. the target is assumed to be intelligent, or at least aware of the interceptor, and tries to maximize the miss distance. In the problem formulation, the pursuer, i.e., interceptor, wishes to minimize the terminal miss whereas the evader, i.e., target, wishes to maximize it. Therefore, the guidance law for the interceptor is determined based on the anticipated worst possible target acceleration, i.e., based on game theoretic principals. Since this guidance law presumes that the engagement states are known by both players, the implementation requires that noisy information be processed through a filter. The filter that is usually chosen also presumes the worst possible target acceleration. However, if the target does not execute according to its assumed strategy, the filter performance can degrade rapidly, and thereby, produce large miss distances. Therefore, there remains a need for a accommodating filtering techniques that may be used with the type of guidance laws based on a game theoretic solution.

SUMMARY OF THE INVENTION

The invention, in its several embodiments of methods and processing structures, generates estimations of intercept and homing missile guidance commands based on the presumption that certain targets execute evasive maneuvers orthogonal to their velocity vectors. The method and processing structure embodiments of the present invention estimate the engagement states in the presence of large and arbitrary target maneuvers and, each processing structure, as a cascade of a filter followed by a guidance law, guide the interceptor, typically by way of acceleration commands orthogonal to the interceptor's velocity vector, in order to hit, or come near to, the target based on these estimated engagement states. One exemplary process and filtering technique constructs the filter by blocking the dynamic effect of the unknown target acceleration and the guidance law based on the anticipated worst possible target acceleration. The relative velocity in the direction of the unknown target acceleration may be constructed directly from the measurements or derived measurements. Another exemplary process and filtering structure constructs the filter that estimates the target acceleration and the guidance law based on a target model in which the target acceleration is assumed to have constant magnitude and constant bank angle rate. The first exemplary process and filtering structure avoids estimating the target acceleration, but must use unfiltered state estimates in the direction of the unknown target acceleration. The second exemplary process and filtering structure works to estimate the target acceleration, but the convergence of this estimate is latent causing miss distance errors.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the present invention in its several embodiments, and for further features and advantages, reference is now made to the following description taken in conjunction with the accompanying drawings, in which:

FIG. 1 is an exemplary top level functional block diagram of the homing missile guidance and estimation processing structure;

FIG. 2 is an exemplary top level functional block diagram of the Target Acceleration Blocking Filter/Game Theoretic Guidance Law (TABF/GTGL) processing structure;

FIG. 3 is an illustration relating an exemplary interceptor position with an exemplary target position;

FIG. 4 is an illustration relating an exemplary interceptor velocity with an exemplary target velocity;

FIG. 5 is an illustration presenting an exemplary target velocity in an inertial coordinate frame;

FIG. 6 is an illustration presenting a target acceleration orthogonal to the target velocity vector;

FIG. 7 is an exemplary top level functional block diagram of the Target Acceleration Estimator Linear Quadratic Guidance Law (TABF/LQGL) processing structure;

FIG. 8 is an exemplary functional block diagram of the filter and guidance law processing structure as part of a fire control of an interceptor; and

FIG. 9 is an exemplary functional block diagram of the filter and guidance law processing structure as part of the guidance, navigation and control system of an interceptor.

DETAILED DESCRIPTION 1. Introduction

The following describes methods and systems for performing the intercept of one vehicle by another. An exemplary embodiment of the system and method embodiments of the present invention may be applied to addresses what is termed in the art as the missile-target intercept problem where one interceptor, or missile, attempts to intercept, or otherwise come into close spatial proximity with, another vehicle such as an incoming weapon that itself acts as an evader. Embodiments of the invention that are the methods and systems have an estimation processing structure and a guidance processing structure that function together to reduce the terminal miss distance between the target and the interceptor.

FIG. 1 illustrates a top level functional block diagram representation of the interceptor system 100 having an interceptor 110 and the ability to account for target dynamics 120. In this example, the interceptor 110 has guidance command and state estimation processing unit 130 that includes the filter processing structure 134 and the guidance law processing structure 132. The filter processing structure 134 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two. The guidance law processing structure 132 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two, or other electrically equivalent means. Accordingly, the filter processing structure 134 and the guidance law processing structure 132, as instructions, may be executed by one or more processors within the guidance command and state estimation processing unit 130. The guidance command and state estimation processing unit 130 takes in measurements 142 from one or more sensors 140 and outputs acceleration commands 138 to the autopilot 160. The autopilot 160 is typically as system within the interceptor 110 that attempts to generate accelerations in the direction and magnitude corresponding to the acceleration commands 138 form the guidance law 132. The autopilot 160 may include control instructions executed via a processor that: (a) takes in angular body rotation measurements and linear acceleration measurements corresponding to the angular velocity and changes in the linear velocity of the interceptor 110; and (b) puts out commands to effectors such as control surface actuators, thrust vector actuators, and/or reaction control jets. The effects of the autopilot and the interceptor's aerodynamics are represented by the block labeled interceptor dynamics 150. The one or more sensors 140 provide the combination of the target dynamics 120 and interceptor dynamics 150 as measurements 142 as well as measurements of interceptor velocity and interceptor acceleration.

1.1 Overview

Two exemplary methods and processing structures for homing missile guidance and estimation based on the assumption that certain targets execute evasive maneuvers orthogonal to their velocity vectors are disclosed. Both exemplary methods and processing structures estimate the engagement states in the presence of large and arbitrary target maneuvers and generate guidance signals, e.g., autopilot acceleration commands, to guide the interceptor to hit the target based on these state estimates.

The first exemplary method and processing structure, the Target Acceleration Blocking Filter (TABF/GTGL), may be derived by applying the linear-exponential-Gaussian (LEG) differential game with different information patterns or its limit, the linear-quadratic-Gaussian (LQG) differential game. In the problem formulation, the interceptor seeks to minimize the terminal miss distance whereas the target seeks to maximize it. Therefore, the guidance law, e.g., the GTGL, for the interceptor is determined based on the anticipated worst possible target motion similar to that of the linear-quadratic pursuit-evasion game (LQPEG). Since LQPEG assumes the engagement states are known by both players, implementation requires that noisy information be processed through a filter. The filter that may be typically applied in the LQPEG assumes the strategy of the target. However, if the target does not play its assumed strategy, then the filter performance can degrade rapidly, and thereby, produce unacceptable miss distances. Instead of assuming the strategy of the target, the filters, e.g., the TABF, in the LEG and LQG differential games block the dynamic effect of the target in the direction of the unknown target acceleration. In effect, this requires the TABF to process its measurements with zero memory in a subspace associated with the direction of the unknown target acceleration. Therefore, the relative velocity in the direction of the unknown target acceleration is constructed directly from the measurements. From the same estimates of the TABF, the guidance law obtained by the game theory described above may now be applied. Since the target acceleration is not estimated, the TABF will not have the estimated lag but will have increased estimation uncertainty because unfiltered state estimates are used in the direction of the unknown target acceleration. Furthermore, since no target model is used in the problem formulation, the filter structure and implementation are much simpler and the computation requirement is a lot less. Finally, since there is no assumption made on the target, the TABF will perform well for a range of different target accelerations. FIG. 2 illustrates in a top level block diagram 200 the TABF/GTGL processing 230 having a target acceleration blocking filter 234 taking in measurements representative of the combined intercept dynamics of the target dynamics 120 and the interceptor dynamics 150, where the exemplary measurements 242 include and elevation angle, el, and azimuth angle, az, relative range, R, and relative range rate, {dot over (R)}. In addition, the interceptor velocity and interceptor acceleration may be provided as measurements. The one or more sensors 140 of the interceptor need not output the exemplary measurements 242, but instead output components, i.e., constituent elements, of the exemplary measurements 242, from which the exemplary measurements 242 may be formed or derived. The game theoretic guidance law 232 processing structure operates on the state estimate 236, {circumflex over (x)}, output by the TABF structure 234, to output the acceleration command 238, u₁, for the autopilot 160.

The TABF processing structure 234 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two, or other electrically equivalent means. The GTGL processing structure 232 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two. Accordingly, the TABF processing structure 234 and the GTGL processing structure 232, as instructions, may be executed by one or more processors within the TABF/GTGL processing unit 230. The TABF/GTGL processing unit 230 takes in measurements 242, in direct or derived form, or combinations thereof, from one or more sensors 140 and outputs acceleration commands 238 to the autopilot 160.

As illustrated by example in FIG. 3, the intercept geometry 300 may be represented by relating the interceptor 310 and the target 320 where an aspect angle 330 may represent the angle between a reference axis 360 and the relative range vector 350, R, and where an azimuth angle 340 may represent the angle between a second reference axis 370, orthogonal to the first reference axis 360, and a projection of the relative range vector 350, R, onto the plane 370 orthogonal to the first reference axis 360. The target 320 may be represented in FIG. 4 as having a velocity vector 410, v_(T), and an acceleration vector, itself having components represented as a_(T) 420 orthogonal to the target velocity vector 410. The interceptor 310 may be represented as having, in a relative frame, a velocity vector that may have a heading error 440 when not aligned with the relative velocity vector projecting an intercept, i.e., the heading error 440 is a deviation from the collision course 450. Illustrated in FIG. 5 is an exemplary inertial frame 500 having three orthogonal axes where the target velocity 510, v_(T), may be represented as a vector v _(T) and expressed in terms of its magnitude 510, v_(T),f light path angle 530, γ, and heading angle 520, Ψ, where the heading angle 520, Ψ, is the angle of rotation from a first axis, e.g., the x-axis, to a projection of the target velocity vector onto the x-z plane, in this example, and where the flight path angle 530, γ, is the angle of rotation from the x-z plane projection to the target velocity vector v _(T).

FIG. 6 illustrates a frame 600 in which the target acceleration 640, α_(T), orthogonal to the target velocity vector 610, V_(T), may expressed in two orthogonal components, e.g., α₁ 620 and a₂ 630 in a plane orthogonal to the velocity vector 610, V_(T). The two orthogonal components may be related to the magnitude of the target acceleration 640, a_(T), by the bank angle 650, φ, where the bank angle 650 is the angular rotation from the axis of the target acceleration component a₁ to the target acceleration vector 640.

The second exemplary method and processing structure, the Target Acceleration Estimator/Linear Quadratic Guidance Law (TAE/LQGL), is derived based on a new target model where the target acceleration vector is assumed to be orthogonal to the target velocity vector and have constant magnitude and constant bank angle rate. This is based on the notion that the target acceleration vector diffuses on a sphere with unknown magnitude, which may be more rational than the usual Gauss-Markov stochastic target model. In addition to representing the actual target maneuvers better, this target model allows a system formulation that results in an implementable filter and guidance law. The extended Kalman filter is used to estimate the engagement states including the magnitude and bank angle of the target acceleration. Alternatively, the modified-gain extended Kalman filter may be used since it allows the measurements of the target in a spherical coordinate system to be processed almost linearly in a rectangular coordinate system. The guidance law is constructed by minimizing a quadratic performance index that includes the terminal miss distance and interceptor maneuverability subject to the assumed engagement dynamics. This guidance law is determined in closed form where the guidance gain is an explicit function of the estimated target acceleration and time-to-go. Furthermore, the LQ guidance law has an adaptive feature in that the guidance gain is adaptive to the target angular rate which is also estimated. Finally, the performance of this second method and processing structure, i.e., the TAE/LQGL, typically depends on whether the assumed target model represents the actual target maneuver as well as the lag in estimating the target acceleration. Note that the first method and processing structure, i.e., the TABF/GTGL, does not typically include estimation lag and may accommodate differing presumed target maneuvers.

FIG. 7 illustrates in a top level block diagram 700 the TAE/LQGL processing 730 having a target acceleration estimator (or filter) 734 taking in measurements representative of the combined intercept dynamics of the target dynamics 120 and the interceptor dynamics 150, where the exemplary measurements 742 include and elevation angle, el, and azimuth angle, az, relative range, R, and relative range rate, {dot over (R)}. In addition, the interceptor velocity and acceleration may be provided as measurements. The one or more sensors 140 of the interceptor need not output the exemplary measurements 742, but instead output components, i.e., constituent elements, of the exemplary measurements 742, from which the exemplary measurements 742 may be formed or derived. The linear quadratic guidance law 732 processing structure operates on the state estimate 736, {circumflex over (x)}, output by the TAE filter structure 234, to output the acceleration command 738, u₁, for the autopilot 160.

The TAE processing structure 734 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two, or other electrically equivalent means. The LQGL processing structure 732 may embodied as instructions executed by a general processor making the general processor function as a special processor or may be embodied as permanent instruction in an application specific integrated circuit or in any combination of the two. Accordingly, the TAE processing structure 734 and the LQGL processing structure 732, as instructions, may be executed by one or more processors within the TAE/LQGL processing unit 730. The TAE/LQGL processing unit 730 takes in measurements 742, in direct or derived form, or combinations thereof, from one or more sensors 140 and outputs acceleration commands 738 to the autopilot 160.

The filter/law processing structures of the TABF/GTGL and TAE/LQGL and the execution of the processes of the TABF/GTGL and TAE/LQGL may be wholly within an interceptor, wholly external to the interceptor, or distributed between the interceptor and a target tracking system for example. That is, the filter/law processing structure may be modular and may be interposed between the integrated measurements and the autopilot whether part of a target tracking system such as a fire control system apart from the interceptor or part of the onboard processing of the interceptor itself. Illustrated in FIG. 8 is an exemplary system in functional block diagram form having a fire control 830, an interceptor 810 and a target 820 wherein the filter/law processing 840 is integrated into the fire control system 830. The target tracking system, a fire control system 830 in this example, includes a target sensor 831 for receiving electromagnetic waves 880 from the target wherein the target sensor 831 may passively sense the target via infrared emissions and/or optical or other radio frequency reflections from the target and the target sensor may actively transmit RF beams to the target. The exemplary fire control system 830 may also sense, via an interceptor sensor 834 or via the target sensor 831, electromagnetic waves 882 from the interceptor 810. Measurements from the target sensor 831 and interceptor sensor 834 may be integrated and/or placed into a reference frame and this may be done via processing steps in a measurement integration module 860. The measurement set 841, e.g., elevation, azimuth, relative range and relative range rate, may then be provided to the filter/law processing module 840. The interceptor velocity and interceptor acceleration may be derived from the interceptor sensor 834 or taken from the inertial measurement sensors or unit that may be integral to an inertial navigation system 813 that may be aided by a Global Positioning Satellite (GPS) receiver system. The interceptor velocity and interceptor acceleration taken from the inertial measurement sensors may be down-linked, typically via an encrypted RF or laser transmission to a downlink receiving system of the fire control 830 and thereafter provided, typically in decrypted form to the filter/law processing module 840. The filter/law processing module 840 may employ a filter/law data store 843, for example, the storing of guidance law gains as a function of anticipated remaining stages in the intercept scenario. Whether the filter law processing module processes the measurements according to the TABG/GTGL processing structure and methods or according to the TAE/LQGL processing structure and methods, the resulting acceleration command 842 may be provided to an interceptor command uplink 832 which may then transmit, typically via an encoded RF or laser link to a command receiver 814 onboard the interceptor 810. The onboard autopilot processing takes in the acceleration command 890, typically decoded by the command receiver 814, as well as measured interceptor dynamics 852 via inertial instruments such as those that may be integral to an inertial navigation system 813 that may be aided via GPS receiver processing. The autopilot generates commands to effectors such as controls surface actuators 812 or reaction control jets or trust vector control actuators. The interceptor dynamics functional block diagram 850 is shown straddling the interceptor 810 to indicate that while not part of the onboard instrumentation and processing, the interceptor dynamics 850 are integral to the intercept scenario.

Illustrated in FIG. 9 is an exemplary system in functional block diagram form having a target tracking system such as an exemplary fire control 930, an interceptor 910 and a target 820 wherein the filter/law processing 940 is integrated into the interceptor 910. The target tracking system, a fire control system 930 in this example, includes a target sensor 931 for receiving electromagnetic waves 980 from the target wherein the target sensor 931 may passively sense the target via infrared emissions and/or optical or other radio frequency reflections from the target and the target sensor may actively transmit RF beams to the target. The exemplary fire control system 930 may transmit target trajectory information including target velocity and acceleration vectors to the interceptor 910 via a direct electrical linked (not shown) if connected and otherwise may transmit target information 984 via a target state uplink 932 to a target state receiver 914 onboard the interceptor 910. The interceptor may have a target sensor 916 may passively sense the target via infrared emissions and/or optical or other radio frequency reflections from the target and the target sensor may actively transmit RF beams to the target. Measurements from the target sensor 931, as may be encoded and uplinked via the target state uplink 932 and received and decoded via the target state receiver 914, and measurements from the interceptor's target sensor 916 may be integrated and/or placed into a reference frame and this may be done via processing steps in a measurement integration module 960 on board the interceptor 910. The measurement set 941, e.g., elevation angle, azimuth angle, relative range and relative range rate, may then be provided to the filter/law processing module 940. The measured velocity and acceleration 951 of the interceptor may be provided to the filter/law processing module 940 via the inertial measurement unit 913. The filter/law processing module 940 may employ a filter/law data store 943, for example, the storing of guidance law gains as a function of anticipated remaining stages in the intercept scenario. Whether the filter law processing module processes the measurements according to the TABG/GTGL processing structure and methods or according to the TAE/LQGL processing structure and methods, the resulting acceleration command 942 may be provided to the onboard autopilot processing 911. The autopilot processing 911 takes in the acceleration command 942, as well as measured interceptor dynamics 952 via inertial instruments such as an inertial measurement unit 913 that may be aided via GPS receiver processing 915. The autopilot generates commands to effectors such as trust vector control actuators 912, controls surface actuators or reaction control jets. The interceptor dynamics functional block diagram 950 is shown straddling the interceptor 910 to indicate that while not part of the onboard instrumentation and processing, the interceptor dynamics 950 are integral to the intercept scenario. An exemplary embodiment where the filter/law processor in onboard the interceptor is used in the simulated results summarized in Tables 1-4.

1.2 Process

Either of the two homing vehicle guidance and estimation processing structures may be described by a high level block diagram when applied against maneuvering targets. There are two major components described as comprising a processing structure and shown in FIG. 1, i.e., a filter that estimates the engagement states and a guidance law that guides the interceptor toward the target. The first exemplary processing structure includes a target acceleration blocking filter (TABF) and a game theoretic guidance law (GTGL). The second exemplary processing structure (TAE/LQGL) includes a target acceleration estimator (TAE) and a linear quadratic guidance law (LQGL).

The inputs to the exemplary process are any combinations of the bearing angles (for example a first bearing angle may be an elevation angle and a second bearing angle may be an azimuth angle), the range and the range rate (or closing velocity) between the target and interceptor measured in an inertia coordinate frame by either a target sensor such as a seeker on the interceptor or any other external sources. The filter combines the inputs through a specialized estimation structure. For the first filtering method and processing structure, the TABF is used. For the second filtering method and processing structure, the TAE is used. The outputs of the filter are typically the estimates of the engagement states. Then, the guidance law combines these state estimates through a specialized guidance structure. For the first guidance method and processing structure, the GTGL is used. For the second guidance method and processing structure, the LQGL is used. The output of the guidance law is typically the interceptor acceleration command which, when effected by the interceptor, will guide the interceptor toward the target. Accordingly, this interceptor acceleration command is the output of the exemplary process which is sent to the autopilot of the interceptor in order to generate the commanded interceptor acceleration.

1.3 Other Applications

While air vehicle interceptors and maneuvering air vehicle targets are employed as exemplary embodiments to teach the invention, the same processes and processing structures are applicable for any type of vehicle attempting to intercept another vehicle.

2. TABF/GTGL Method Processing and Processing Structure

For purposes of describing the TABF/GTGL processing, one may consider the three-dimensional relative dynamics between the target and interceptor in an inertial coordinate frame. States may be defined in the following set of equations as:

-   -   {dot over (x)}_(r)=u_(r)     -   {dot over (y)}_(r)=u_(r)     -   ż_(r)=ω_(r)         {dot over (u)} _(r) =a _(1x) −a _(1x)  [1]         {dot over (v)} _(r) =a _(1y) −a _(1y)         {dot over (ω)}_(r) =a _(Tz) −a _(1z)         where x, y, and z, represent components of the relative position         vector, u_(r), v_(r) and ω_(r) represent components of the         relative velocity vector, a_(Tx), a_(Ty) and a_(tz) represent         components of the target acceleration vector and a_(1x), a _(1y)         and a_(1z) represent components of the interceptor acceleration         vector. There are at least four possible measurements, where         four are used in this example for the interceptor having one or         more instruments that measure directly, or determines, the         elevation angle, el, the azimuth angle, az, the range, R and the         range rate, {dot over (R)}, between the target and interceptor         and they be related in the following set of equations:

$\begin{matrix} \begin{matrix} {{el} = {\sin^{- 1}\frac{y_{r}}{\sqrt{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}}}} \\ {{az} = {\tan^{- 1}\frac{z_{r}}{x_{r}}}} \\ {R = \sqrt{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}} \\ {\overset{.}{R} = \frac{{x_{r}u_{r}} + {y_{r}\upsilon_{r}} + {z_{r}\omega_{r}}}{\sqrt{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}}} \end{matrix} & \lbrack 2\rbrack \end{matrix}$

Note that here the elevation angle is defined with respect to the xz-plane and the azimuth angle is defined in the xz-plane. However, these two angles can also be defined using the xy-plane or yz-plane. In addition, the interceptor has measurements of its own velocity, denoted as v_(1x), v_(1y) and v_(1z), and acceleration, i.e., a_(1x), a_(1y) and a_(1z).

An objective of the missile intercept problem is to determine the interceptor acceleration such that the terminal miss distance is minimized subject to the relative dynamics of the set of equation set 1 taking in all four measurements in the equation set 2 or any subset of the four measurements. This problem is solved by applying the solution of the linear-exponential-Gaussian (LEG) differential game or its limit, the linear-quadratic-Gaussian (LQG) differential game. In Section 2.1, the system dynamics supporting the TABF/GTGL processing structure are described. In Section 2.2, exemplary solutions of the LEG and LQG differential games are described. In Section 2.3, the generalizations for applying the solutions of the LEG and LQG differential games to the air vehicle, particularly missile target, intercept problem are described. In Section 2.4, the TABF/GTGL processing structure is summarized. In Section 2.5, an illustrative numerical example is provided.

2.1 System Dynamics

State and control variables may be represented as follows:

$\begin{matrix} {{x = \begin{bmatrix} x_{r} \\ y_{r} \\ z_{r} \\ u_{r} \\ \upsilon_{r} \\ \omega_{r} \end{bmatrix}},{u_{T} = \begin{bmatrix} a_{Tx} \\ a_{Ty} \\ a_{Tz} \end{bmatrix}},{u_{I} = \begin{bmatrix} a_{Ix} \\ a_{Iy} \\ a_{Iz} \end{bmatrix}}} & \lbrack 3\rbrack \end{matrix}$

Then, equation set 1 with the addition of process noise ω can be written as {dot over (x)}=Ax+B ₁ u _(T) +B ₁ u ₁+ω  [4] where

$\begin{matrix} {{A = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}},{B_{T} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}},{B_{i} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ {- 1} & 0 & 0 \\ 0 & {- 1} & 0 \\ 0 & 0 & {- 1} \end{bmatrix}}} & (5) \end{matrix}$

Continuous time dynamics may be represented by the propagations of finite time increments, or discrete time intervals, via a discretization of the continuous time dynamics. By discretizing equation 4, the discrete-time relative dynamics may be represented as: x _(i+1) =φx _(i) +Γ ₁ u _(T,i) +Γ _(l) u _(l.i) +ω _(i)  [6] where

$\begin{matrix} {{\Phi = {I + {A\;\Delta\; t}}}{\Gamma_{T} = {\left( {I + \frac{A\;\Delta\; t}{2}} \right)B_{T}\Delta\; t}}{\Gamma_{I} = {\left( {I + \frac{A\;\Delta\; t}{2}} \right)B_{I}\Delta\; t}}} & \lbrack 7\rbrack \end{matrix}$ and where Δt is the sampling time. From the system of equation 2, the measurement equation with the addition of sensor noise v_(i) can be written as z _(i) =h(x _(i))+v _(i)  [8] where z_(i) may include four measurements, i.e., two angles, range and range rate, or an subset of the four measurements. It is assumed for these exemplary embodiments that ω_(i) and u_(i) are zero mean, white Gaussian process with variance W and V, respectively. In addition, it assumed for these exemplary embodiments that the interceptor has measurements of its own velocity, i.e., v_(1x), v_(1y) and v_(1z), and acceleration, i.e., a_(1x), a_(1y) and a_(1z), available to its processing.

2.2 Solutions of LEG and LQG Differential Games

The LEG differential game may be expressed as:

$\begin{matrix} {\min\limits_{u_{I}}{\max\limits_{u_{T}}{E\left\lbrack e^{{- \gamma}\; S} \right\rbrack}}} & \lbrack 9\rbrack \end{matrix}$ where

$\begin{matrix} {S = {{\frac{1}{2}{{x_{N}}}_{Q_{N}}^{2}} + {\sum\limits_{i = 0}^{N - 1}\;{\frac{1}{2}\left( {{{{u_{I,i}}}_{R_{I}}^{2} +}❘{u_{T,i}{}_{R_{T}}^{2}}} \right)}}}} & \lbrack 10\rbrack \end{matrix}$ subject to x _(i+1) =Φx _(i) +Γ _(T) u _(T,i) +Γ ₁ u _(1,i)+ω_(i) z _(i) =H _(i) x _(i) +v _(i)  [11] where N is the number of stages known a priori, y≧0, Q_(N)≧0, R₁>0 and R_(T)<0 are design parameters. As γ→0, the LEG differential game becomes the LQG differential game. The solution of the LEG differential game, as a processing structure, includes a filter and a guidance law. The filter includes the time propagation and measurement update steps. The time propagation of the state and covariance may be represented as: x _(i+1) =Φ{circumflex over (x)} _(i) +Γu _(1,i) P _(i+1) =ΦP _(i)Φ^(T) +W−Γ _(T)(γR _(T))⁻¹Γ_(T) ^(T)  [12] where x is the a prior state estimate, {circumflex over (x)} is the a posteriori state estimate, P is the a prior error covariance and P is the a posteriori error covariance. The measurement update may be processed as {circumflex over (x)} _(i) = x _(i) +K(z _(i) −H _(i) x _(i)) P _(i)=(I−K _(i) H _(i)) P (I−K _(i) H _(i))^(T) +K _(i) VK _(T) ^(i)  [13] And where the filter gain may be processed as K _(i) = P _(i) H _(i) ^(T)(H^(i) P _(i) H _(i) ^(T) +V) ⁻¹.  [14]

Since this filter cannot be evaluated in the limit where γ→0, the time propagation and measurement update are rewritten in terms of P ⁻¹ instead of P. The time propagation may be rewritten and processed as: x _(i+1) =Φ{circumflex over (x)} _(i)+Γ₁ u _(1,i) P _(i+1) ⁻¹ =N _(i) ⁻¹ −N _(i) ⁻¹Φ(P _(i) ⁻¹ +Φ ^(T) N _(i) ⁻¹Φ)Φ^(T) N _(i) ⁻¹  [15] N _(i) ⁻¹ =W ⁻¹ −W ⁻¹Γ_(T)(Γ_(T) ^(T) W ⁻¹ Γ _(T) −γR _(T))⁻¹ Γ _(T) ^(T) W ⁻¹ or x _(x+1) =Φ{circumflex over (x)} _(i) +Γ _(i) u _(1,i) P _(i+1) ⁻¹ =N _(i) ³¹ ¹ −N _(i) ⁻¹Γ_(T)(Γ_(T) ^(T) N _(i) ⁻¹ Γ _(T) −γR _(T))⁻¹ Γ _(T) ^(T) N _(i) ⁻¹  [16] N _(i) ⁻¹=(ΦP ₁ Φ ^(T) +W)⁻¹ The measurement update may be reformulated and processed as {circumflex over (x)} _(i) = x+P _(i) H _(i) ^(T) V ⁻¹(z _(i) −H _(i) x _(i)) P _(i)=( P _(i) ⁻¹ +H _(i) ^(T) V ⁻¹ H _(i))⁻¹  [17]

As previously described, in the limit, i.e., when γ→0, the LEG differential game becomes the LQG differential game. Therefore, the time propagation for the LQG differential game may be processed as: x _(i+1) =Φ{circumflex over (x)} _(i)+Γ₁ u _(1,i) P _(i+1) ⁻¹ =N _(i) ⁻¹ −N _(i) ⁻¹Φ(P _(i) ⁻¹ +Φ ^(T) N _(i) ⁻¹Φ)Φ^(T) N _(i) ⁻¹  [18] N _(i) ⁻¹ =W ⁻¹ −W ⁻¹Γ_(T)(Γ_(T) ^(T) W ⁻¹Γ_(T))⁻¹ Γ _(T) ^(T) W ⁻¹ or x _(i+1) =Φ{circumflex over (x)} _(i) +Γ ₁ u _(1,i) P _(i+1) ⁻¹ =N _(i) ⁻¹ −N _(i) ⁻¹Γ_(T)(Γ_(T) ^(T) N _(i) ⁻¹Γ_(T))⁻¹ Γ _(T) ^(T) N _(i) ⁻¹  [19] N _(i) ⁻¹=(ΦP _(i)Φ^(T) +W)⁻¹

The measurement update for the LQG differential game may be processed as {circumflex over (x)} _(i) = x _(i) +P _(i) H _(i) ^(T) V ⁻¹(z _(i) −H _(i) x _(i)) P _(i)=( P _(i) ⁻¹ +H _(i) ^(T) V ⁻¹ H _(i))⁻¹  [20]

For both LEG and LQG differential games, the guidance law may be expressed as: u _(1,i)=G_(i) {circumflex over (x)} ₁;  [21] where the guidance gain may be expressed as: G _(i) =−R ₁ ⁻¹Γ₁ ^(T)Π_(i+1)(I+Γ ₁ R ₁ ⁻¹Γ₁ ^(T)Π_(i+1) +Γ _(T) R _(T) ⁻¹Γ_(T) ^(T)Π_(i+1))⁻¹Φ;  [22] and where Π_(i)=Φ^(T)Π_(i+1)(I+Γ ₁ R ₁ ⁻¹Γ₁ ^(T)Π_(i+1) +Γ _(T) R _(T) ⁻¹Γ_(T) ^(T)Π_(i+1))⁻¹Φ, Π_(N) =Q _(N)  [23]

2.3 Generalization for the Missile Intercept Problem

In order to apply the solutions of the LEG AND LQG differential games to the missile intercept problem, three generalizations are required. First, the LEG and LQG differential games assume that the number of stages is fixed and known a priori while the number of stages is unknown in the missile intercept problem. Therefore, a generalization may be made to determine the number of remaining stages to be used in the guidance law. Since the time-to-go may be approximated by

$\begin{matrix} {T_{go} = {\frac{R}{\overset{.}{R}} = \frac{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}{{x_{r}u_{r}} + {y_{r}\upsilon_{r}} + {z_{r}\omega_{r}}}}} & \lbrack 24\rbrack \end{matrix}$ where x_(r), y_(r), z_(r), u_(r), v_(r) and ω_(r) are estimated as {circumflex over (x)}_(i), the number of remaining stages can be estimated online by dividing the estimated times-to-go by the sampling time Δt and rounding to the nearest integer.

Second, the LEG and LQG differential games typically apply linear measurements while the measurements of equation 8 used in the missile intercept problem are nonlinear. Therefore, a generalization may be required to modify the filter to a structure similar to that of an extended Kalman filter. While the time propagation remains the same, the measurement update for both LEG and LQG differential games in the exemplary scenario becomes {circumflex over (x)}= x _(i) +K _(i) [z _(i) −h( x _(i))] P _(i)=(I−K _(i) H _(i)) P _(i)(I−K _(i) H _(i))^(T) +K _(i) VK _(i) ^(T)  [25] K _(i) = P _(i) H _(i) ^(T)(H _(i) P _(i) H _(i) ^(T) +V)⁻¹ or {circumflex over (x)} _(i) = x _(i) +P _(i) H _(i) ^(T) V ⁻¹ [z _(i) −h( x _(i))]  [26] P _(i)=( P _(i) ⁻¹ +H _(i) ^(T) V ⁻¹ H _(i))⁻¹ where

$\begin{matrix} {H_{i} = {\frac{\partial h}{\partial x_{i}}❘_{x_{i} = {\overset{\_}{x}}_{i}}\mspace{14mu}{and}}} & \lbrack 27\rbrack \\ {{\frac{\partial{el}}{\partial x_{i}} = \left\lbrack {\frac{{- x_{r}}y_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)\left( {x_{r}^{2} + z_{r}^{2}} \right)^{0.5}}\frac{\left( {x_{r}^{2} + z_{r}^{2}} \right)^{0.5}}{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}\mspace{14mu}\ldots\mspace{14mu}\frac{{- y_{r}}z_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)\left( {x_{r}^{2} + z_{r}^{2}} \right)^{0.5}}000} \right\rbrack}{\frac{\partial{az}}{\partial x_{i}} = \left\lbrack {\frac{- z_{r}}{x_{r}^{2} + z_{r}^{2}}0\frac{x_{r}}{x_{r}^{2} + z_{r}^{2}}000} \right\rbrack}{\frac{\partial R}{\partial x_{i}} = \left\lbrack {\frac{x_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}\frac{y_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}\frac{z_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}000} \right\rbrack}{\frac{\partial\overset{.}{R}}{\partial x_{i}} = \left\lbrack {\frac{{u_{r}\left( {y_{r}^{2} + z_{r}^{2}} \right)} - {x_{r}\left( {{y_{r}\upsilon_{r}} + {z_{r}\omega_{r}}} \right)}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{1.5}}\frac{{\upsilon_{r}\left( {x_{r}^{2} + z_{r}^{2}} \right)} - {y_{r}\left( {{x_{r}u_{r}} + {z_{r}\omega_{r}}} \right)}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{1.5}}\mspace{14mu}\ldots\mspace{14mu}\frac{{\omega_{r}\left( {x_{r}^{2} + y_{r}^{2}} \right)} - {z_{r}\left( {{x_{r}u_{r}} + {y_{r}\upsilon_{r}}} \right)}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{1.5}}\frac{x_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}\frac{y_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}\mspace{14mu}\ldots\mspace{14mu}\frac{z_{r}}{\left( {x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} \right)^{0.5}}} \right\rbrack}} & \lbrack 28\rbrack \end{matrix}$

Third, the LEG and LQG differential games typically consider a three-dimensional direction Γ_(T) for the target acceleration while the missile intercept problem may use a two-dimensional direction that is orthogonal to the target velocity vector. Since the target acceleration vector may be expressed as being orthogonal to the target velocity vector, a _(Tx) v _(Tx) +a _(Ty) v _(Ty) +a _(Tz) v _(Tz)=0  [29] where v_(Tx), v_(Ty) and v_(Tz) are the target velocity. Then, a _(Tz) =−v _(Tx) v _(Tz) ⁻¹ a _(Tx) −v _(Ty) v _(Tz) ⁻¹ a _(Ty)  [30] and

$\begin{matrix} {{B_{T}u_{T}} = {{\begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{Tx} \\ a_{Ty} \\ a_{Tz} \end{bmatrix}} = {\begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ \upsilon_{Tz} & 0 \\ 0 & \upsilon_{Tz} \\ {- \upsilon_{Tx}} & {- \upsilon_{Ty}} \end{bmatrix}\begin{bmatrix} {\upsilon_{Tz}^{- 1}a_{Tx}} \\ {\upsilon_{Tz}^{- 1}a_{Ty}} \end{bmatrix}}}} & \lbrack 31\rbrack \end{matrix}$

Now the new two-dimensional target acceleration direction can be estimated as

$\begin{matrix} {\Gamma_{T} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ {\omega_{r} + \upsilon_{Iz}} & 0 \\ 0 & {\omega_{r} + \upsilon_{Iz}} \\ {- \left( {u_{r} + \upsilon_{Ix}} \right)} & {- \left( {\upsilon_{r} + \upsilon_{Iy}} \right)} \end{bmatrix}} & \lbrack 32\rbrack \end{matrix}$ where u_(r), v_(r) and ω_(r) are estimated as part of {circumflex over (x)}_(i) and u_(1x), u_(1y) and u_(1z) are measured and known by the interceptor. Therefore, the target acceleration direction Γ_(T) can be estimated online, i.e., by the processing structure typically in real time, as two-dimensional and time-varying instead of three-dimensional, time-invariant and a priori known. This new target acceleration direction will be used in the time propagation of the filter while the measurement update of the filter does not require the use of target acceleration direction. Furthermore, because this new target acceleration direction is not known a priori, in the revised processing structure, it is not incorporated into the structure of the exemplary guidance law.

2.4 Summary of Generating a Filter/Law Processing Structure

The structuring and implementation procedures of the TABF/GTGL processing structure shown in FIG. 2 are summarized below.

2.4.1 Structuring Procedure

The structuring procedure for the GTGL may be as follows:

-   -   1. Define a sampling time, Δt, which may be derived from a trade         of computational requirements and acceptable propagation error.     -   2. Generate discrete-time dynamics, i.e., Φ, Γ_(T) and Γ₁.     -   3. Choose an arbitrarily large value for N that is sufficient in         value to accommodate the total iterations or stages, expected in         the execution of an intercept.     -   4. Determine structure parameters R₁, R_(T) and Q_(N)         experimentally, that is, to choose these parameters according to         the acceleration capability of the interceptor, the expected         acceleration capability of the target and the desired terminal         miss distance, and to tune these parameters in numerical         simulation.     -   5. Calculate guidance gain as in equation 22 and store it as a         table as a function of the number of the remaining stages.

The structuring procedure for the TABF may be as follows:

-   -   1. Define a sampling time, Δt, which may be derived from a trade         of computational requirements and acceptable propagation error.     -   2. Generate discrete-time dynamics, i.e., Φ and Γ₁.     -   3. Configure the structure to take in at least one of four         measurements from two angles, range and range rate.     -   4. Determine structural parameters W, V, γ and R_(T)         experimentally, that is, to choose these parameters according to         the measurement uncertainty and system disturbances, and to tune         these parameters in numerical simulation.     -   5. Determine an initial P experimentally from the observed state         uncertainties and to account for particular dynamical         uncertainties as process noise.

Note that R_(T) for the filter is different from the R_(T) for the guidance law because the Γ_(T) for the filter is two-dimensional, time-varying and estimated online while the Γ_(T) for the guidance law is three-dimensional, time-invariant and known a priori. Also note that the filter gain may be determined online, that is, by the executing of the processing structure, because of the linearization of h around x and the estimation of Γ_(T) by using {circumflex over (x)}.

2.4.2 Implementation Procedure

At every sampling time, the following steps may be performed:

-   -   1. Acquire the input of the TABF/GTGL processing structure from         the sensors chosen in the design procedure and the sensors that         measure the interceptor's own velocity and acceleration.     -   2. Perform filter measurement update.         -   Determine a linearized measurement matrix as in equation 27             by using the a priori state estimate.         -   Determine the a posteriori state estimate and the a             posteriori error covariance as in equation 25 or equation 26             by using the a priori state estimate, the a priori error             covariance and the acquired measurements.     -   3. Generate interceptor acceleration command.         -   Determine a time-to-go estimate as in equation 24 by using             the a posteriori state estimate.         -   Determine the number of remaining stages by applying the             estimated time-to-go value.         -   Obtain guidance gain through table lookup by using the             number of remaining stages.         -   Determine interceptor acceleration command as in equation 21             by applying the guidance gain and the a posteriori state             estimate.     -   4. Send interceptor acceleration command as the output of the         TABF/GTGL processing structure to the autopilot of the         interceptor to generate interceptor acceleration, e.g., by         re-orienting the interceptor to create a change in the         angle-of-attack of the interceptor that produces the commanded         interceptor acceleration.     -   5. Perform filter time propagation.         -   Determine Γ_(T) as in equation 32 by using the a posteriori             state estimate and interceptor velocity measurement.         -   Determine the a priori state estimate and the a priori error             covariance by using the a posteriori state estimate, the a             posteriori error covariance and interceptor acceleration             measurement as in equation 12, equation 15, or equation 16,             if the LEG differential game solution is used or equation 18             or equation 19 if the LQG differential game solution is             used.

2.5 Numerical Example

In this section, a numerical example is given to demonstrate the performance of the TABF/GTGL method and processing structure. This example simulates a ballistic air vehicle target in an intercept scenario as shown in FIG. 3 and FIG. 4. The initial distance between the target and interceptor is 50 km. The initial aspect angle is 0, 15 or 30 deg. The initial azimuth is 0, 45, 90, 135, 180, 225, 270 or 315 deg. The initial target velocity is 3 km/s pointing toward the ground. The initial interceptor velocity is 2 km/s with the initial heading angle being 0, 5 or 10 deg off from the collision course associated with a non-maneuvering target. The target acceleration vector is orthogonal to the target velocity vector. There are three different target acceleration vectors starting at range-to-go of 20 km. All three target acceleration factors have a magnitude of 7 g, i.e., seven time the acceleration of gravity, while one does not rotate in the plane that is orthogonal to the target velocity vector and the other two rotate at 0.1 and 2 Hz, respectively. The magnitude of the interceptor acceleration vector is limited to be below 20 g and has a lag in its response to changes in acceleration command represented in this example as a first-order filter with time constant of 0.1 sec. There are three different measurement sets. The first set only has angle measurement. The second set has angle and range measurements. The third set has angle, range and range rate measurements. The measurements are updated at 100 Hz. There are two levels of sensor noise in this example. The larger sensor noise has standard deviation of 0.001 rad for the angle measurement, 1 m for the range measurement and 1 m/s for the range rate measurement. The smaller sensor noise has standard deviation of 0.0001 rad for the angle measurement, 0.1 m for the range measurement and 0.1 m/s for the range rate measurement.

The exemplary TABF/GTGL processing is structured for the three measurement sets with two levels of sensor noise and evaluated under the three target accelerations using Monte-Carlo analysis. For each case, 100 simulation runs are conducted and 90% of the data are used to calculate the averaged miss distance. The averaged miss distance of the algorithm when the larger sensor noise is used is summarized in Table 1. The averaged miss distance of the algorithm when the smaller sensor noise is used is summarized in Table 2. In both Tables 1 and 2, the first, second and third rows represent the three cases where the target acceleration does not rotate, rotates at 0.1 Hz and rotates at 2 Hz, respectively. The first, second and third columns represent the three cases where different measurement sets are used. These statistical miss distance results show that the TABF/GTGL processing structure has hit-to-kill capability for all three target accelerations under the ballistic missile intercept scenario.

TABLE 1 Simulated averaged miss distance for interceptor having lower quality sensor modeled. Measurement Angles, Range & Target maneuver Angles Angles & Range Range Rate No rotation 1.3501 m 1.2858 m 0.7182 m Rotate at 0.1 Hz 1.3769 m 1.3091 m 0.7822 m Rotate at 2 Hz 0.6653 m 0.6254 m 0.4630 m

TABLE 2 Simulated averaged miss distance for interceptor having higher quality sensor modeled. Measurement Angles, Range & Target maneuver Angles Angles & Range Range Rate No rotation 0.0940 m 0.0851 m 0.0581 m Rotate at 0.1 Hz 0.0944 m 0.0894 m 0.0591 m Rotate at 2 Hz 0.1007 m 0.0929 m 0.0628 m

3. TAE/LQGL Method and Processing Structure

The following describes the three-dimensional relative dynamics between the target and interceptor in an inertial coordinate frame: {dot over (x)}_(r)=u_(r)  [33.1] {dot over (y)}_(r)=v_(r)  [33.2] ż_(r)=ω_(r)  [33.3] {dot over (u)} _(r) =a _(1x) −ā _(1x)  [33.4] {dot over (u)} _(r) =a _(1x) −ā _(1y)  [33.5] {dot over (ω)} _(r) =a _(1z) −ā _(1z)  [33.6] {dot over (a)} _(1x) =−aā _(1x) +aa _(1x)  [33.7] {dot over (a)} _(1y) =−aā _(1y) +aa _(1y)  [33.8] {dot over (a)} _(1z) =−aā _(1z) +aa _(1z)  [33.9] where x_(r), y_(r) and z_(r) are the relative position, u_(r), v_(r) and ω_(r) are the relative velocity, a_(tx), a_(ty) and a_(1z) are the target acceleration, ā_(1x), ā_(1y) and ā_(1z) are the interceptor acceleration, a_(1x), a_(1y), and a_(1z) are the interceptor acceleration command and a⁻¹ is the known time constant of the interceptor's actuator lag. There are four possible measurement for the interceptor that measure the elevation angle el, azimuth angle az, range R and range rate {dot over (R)} between the target and interceptor, and they may be represented as follows:

$\begin{matrix} {{z_{el} = {{\sin^{- 1}\frac{y_{r}}{\sqrt{x_{5}^{2} + y_{r}^{2} + z_{r}^{2}}}} + \upsilon_{el}}}{z_{az} = {{\tan^{- 1}\frac{z_{r}}{x_{r}}} + \upsilon_{az}}}{z_{R} = {\sqrt{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}} + \upsilon_{R}}}{z_{\overset{.}{R}} = {\frac{{x_{r}u_{r}} + {y_{r}\upsilon_{r}} + {z_{r}\omega_{r}}}{\sqrt{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}} + \upsilon_{\overset{.}{R}}}}} & \lbrack 34\rbrack \end{matrix}$ where v_(el), v_(az), v_(R) and v_({dot over (R)}) are zero mean, white Gaussian processes with variances V_(el), V_(az), V_(R) and V_({dot over (R)}), respectively, representing the measurement uncertainty. Note that here the elevation angle is defined with respect to the xz-plane and the azimuth angle is defined in the xz-plane. However, these two angles can also be defined using xy-plane or yz-plane. In addition, the interceptor has measurements of its own velocity, denoted as u_(1x), u_(1y) and u_(1z), and acceleration. i.e., ā_(1x), ā_(1y) and ā_(1z).

The objective of the missile intercept problem is typically to determine the interceptor acceleration command, or an equivalent, such that the terminal miss distance is minimized subject to the relative dynamics of equations 33.1-33.9 given all four measurements in the set of equations 34 or any subset of the four measurements. This problem may be solved by using a target acceleration estimator (TAE) and a linear quadratic guidance law (LQGL). It may be assumed in the present exemplary solution that the target acceleration vector s orthogonal to the target velocity vector. Furthermore, it may be assumed that the magnitude and bank angle rate of the target acceleration vector are constant. The TAE estimates the magnitude and bank angle of the target acceleration vector based on a target model and the LQGL is designed by minimizing a quadratic performance index that includes terminal miss distance and interceptor maneuverability subject to the assumed engagement dynamics.

3.1 Target Model

In this section, the target model, or a representation of the target's dynamical characteristics, is developed based on the exemplary assumptions that the target acceleration vector s orthogonal to the target velocity vector and the magnitude and bank angle rate of the target acceleration vector are constant. First, the target velocity vector v _(T) is expressed in terms of its magnitude v_(T), flight path angle γ and heading angle Ψ as v _(T) =v _(T) cos γ cos Φ i +v _(T) sin γ j+v _(T) cos γ sin Ψ k   [35] where ī, j and k are the unit vectors along the three axis of an inertial coordinate frame as in FIG. 5. In order to define the plane for the target acceleration vector ā_(T) that is orthogonal to the target velocity vector, two unit vectors are needed. One vector can be obtained by rotating the unit vector along v _(T) by 90 degrees in γ as ā ₁=−sin γ cos Ψī+cos γ j −sin γ sin Ψ k.  [36]

The other vector can be obtained as ā ₂ =ā ₁ ×ū _(T) =sin Ψī−cos Ψ k.  [37]

Therefore, ā_(T) can be expressed as ā _(T) =a _(T) cos φā ₁ +a _(t) sin φ a ₂,  [38] where a_(T) and φ are the magnitude and bank angle of the target acceleration vector, respectively, as shown in FIG. 6.

By substituting equation 36 and equation 37 into equation 38, ā _(T)=(−a _(T) cos φ sin γcos Ψ+a _(T) sin φ sin Ψ)ī+a _(T) cos φ cos γ j+(−a _(T) cos φ sin γ sin Ψ−a _(T) sin φ cos Ψ) k   [39 ]

By using Brownian motion processes ω_(γ) and ω_(ψ) to represent the uncertainty in the orthogonality between the target acceleration and velocity vectors, equation 39 becomes ā _(T) =[−a _(T) cos φ sin(γ+ω_(γ)) cos (Ψ+ω_(Ψ))+a _(T) sin φ sin (Ψ+ω_(Ψ))]ī+a _(T) cos φ cos (γ+ω_(γ)) j+[−a _(T) cos φ sin (γ+ω_(γ)) sin (ψ+ω_(ψ))−a _(T) sin φ cos (ψ+ω_(ψ))] k   [40]

By defining the target states as: a ₁ =a _(T) sin φ sin (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₂ =a _(T) sin φ sin (γ+ω_(γ)) cos (ψ+ω_(ψ)) a ₃ =a _(T) sin φ cos (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₄ =a _(T) sin φ cos (γ+ω_(γ)) cos (ψ+ω_(ψ)) a ₅ =a _(T) cos φ sin (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₆ =a _(T) cos φ sin (γ+ω_(γ)) cos (ψ+ω_(ψ)) a ₇ =a _(T) cos φ cos (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₈ =a _(T) cos φ cos (γ+ω_(γ)) cos (ψ+ω_(ψ)) a ₉ =a _(T) sin φ sin (γ+ω_(γ)) a ₁₀=a_(T) sin φ cos (γ+ω_(γ)) a ₁₁ =a _(T) cos φ sin (γ+ω_(γ)) a ₁₂ =a _(T) cos φ cos (γ+ω_(γ)) a ₁₃ =a _(T) sin φ sin (ψ+ω_(ψ)) a ₁₄ =a _(T) sin φ cos (ψ+ω_(ψ)) a ₁₅ =a _(T) cos φ sin (ψ+ω_(ψ))  [41] a ₁₆ =a _(T) cos φ cos (ψ+ω_(ψ)) then, equation 40 becomes ā _(T)=(−a ₆ +a ₁₃)ī+a ₁₂ j+(−a ₅ −a ₁₄) k   [42] and accordingly a _(1x) =−a ₆ +a ₁₃  [43] a_(1y=a) ₁₂ a _(1z) =−a ₅ −a ₁₄ may be used to express the target acceleration vector components.

By taking the Ito stochastic differential of the system of equation 41, the target model is obtained as

$\begin{matrix} {\begin{bmatrix} {\overset{.}{a}}_{1} \\ {\overset{.}{a}}_{2} \\ {\overset{.}{a}}_{3} \\ {\overset{.}{a}}_{4} \\ {\overset{.}{a}}_{5} \\ {\overset{.}{a}}_{6} \\ {\overset{.}{a}}_{7} \\ {\overset{.}{a}}_{8} \end{bmatrix} = {{\begin{bmatrix} {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \overset{.}{\psi} & \overset{.}{\gamma} & 0 & \overset{.}{\phi} & \ldots & 0 & 0 & 0 \\ {- \overset{.}{\psi}} & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & 0 & \overset{.}{\gamma} & 0 & \ldots & \overset{.}{\phi} & 0 & 0 \\ {- \overset{.}{\gamma}} & 0 & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \overset{.}{\psi} & 0 & \ldots & 0 & \overset{.}{\phi} & 0 \\ 0 & {- \overset{.}{\gamma}} & {- \overset{.}{\psi}} & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & 0 & \ldots & 0 & 0 & \overset{.}{\phi} \\ {- \overset{.}{\phi}} & 0 & 0 & 0 & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \ldots & \overset{.}{\psi} & \overset{.}{\gamma} & 0 \\ 0 & {- \overset{.}{\phi}} & 0 & 0 & {- \overset{.}{\psi}} & \ldots & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & 0 & \overset{.}{\gamma} \\ 0 & 0 & {- \overset{.}{\phi}} & 0 & {- \overset{.}{\gamma}} & \ldots & 0 & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \overset{.}{\psi} \\ 0 & 0 & 0 & {- \overset{.}{\phi}} & 0 & \ldots & {- \overset{.}{\gamma}} & {- \overset{.}{\psi}} & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} \end{bmatrix}\begin{bmatrix} \begin{matrix} a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \\ a_{5} \\ a_{6} \\ a_{7} \end{matrix} \\ a_{9} \end{bmatrix}} + {\begin{bmatrix} a_{3} & a_{2} \\ a_{4} & {- a_{1}} \\ {- a_{1}} & a_{4} \\ {- a_{2}} & {- a_{3}} \\ a_{7} & a_{6} \\ a_{8} & {- a_{5}} \\ {- a_{5}} & a_{8} \\ {- a_{6}} & {- a_{7}} \end{bmatrix}\begin{bmatrix} {\overset{.}{\omega}}_{\gamma} \\ {\overset{.}{\omega}}_{\psi} \end{bmatrix}}}} & \lbrack 44.1\rbrack \\ {\begin{bmatrix} {\overset{.}{a}}_{9} \\ {\overset{.}{a}}_{10} \\ {\overset{.}{a}}_{11} \\ {\overset{.}{a}}_{12} \end{bmatrix} = {{\begin{bmatrix} {- \frac{\Omega_{\gamma}}{2}} & \overset{.}{\gamma} & \overset{.}{\phi} & 0 \\ {- \overset{.}{\gamma}} & {- \frac{\Omega_{\gamma}}{2}} & 0 & \overset{.}{\phi} \\ {- \overset{.}{\phi}} & 0 & {- \frac{\Omega_{\gamma}}{2}} & \overset{.}{\gamma} \\ 0 & {- \overset{.}{\phi}} & {- \overset{.}{\gamma}} & {- \frac{\Omega_{\gamma}}{2}} \end{bmatrix}\begin{bmatrix} a_{9} \\ a_{10} \\ a_{11} \\ a_{12} \end{bmatrix}} + {\begin{bmatrix} a_{10} \\ {- a_{9}} \\ a_{12} \\ {- a_{11}} \end{bmatrix}{\overset{.}{\omega}}_{\gamma}}}} & \lbrack 44.2\rbrack \\ {{{\begin{bmatrix} {\overset{.}{a}}_{13} \\ {\overset{.}{a}}_{14} \\ {\overset{.}{a}}_{15} \\ {\overset{.}{a}}_{16} \end{bmatrix}\begin{bmatrix} {- \frac{\Omega_{\psi}}{2}} & \overset{.}{\psi} & \overset{.}{\phi} & 0 \\ {- \overset{.}{\psi}} & {- \frac{\Omega_{\psi}}{2}} & 0 & \overset{.}{\phi} \\ {- \overset{.}{\phi}} & 0 & {- \frac{\Omega_{\psi}}{2}} & \overset{.}{\psi} \\ 0 & {- \overset{.}{\phi}} & {- \overset{.}{\psi}} & {- \frac{\Omega_{\psi}}{2}} \end{bmatrix}}\begin{bmatrix} {\overset{.}{a}}_{13} \\ {\overset{.}{a}}_{14} \\ {\overset{.}{a}}_{15} \\ {\overset{.}{a}}_{16} \end{bmatrix}} + {\begin{bmatrix} a_{14} \\ {- a_{13}} \\ a_{16} \\ {- a_{15}} \end{bmatrix}{\overset{.}{\omega}}_{\psi}}} & \lbrack 44.3\rbrack \end{matrix}$ with {umlaut over (φ)}={dot over (ω)}_(φ)  [44.4] and where {dot over (ω)}_(γ) and {dot over (ω)}_(ψ) are zero mean, white Gaussian processes with power spectral densities Ω_(γ) and Ω_(ψ), respectively, and {dot over (ω)}_(φ) is a zero mean, white Gaussian process with power spectral density Ω_(φ) representing the uncertainty in the constant bank angle rate. By taking the derivative of equation 35 and comparing to equation 39, the time derivative of the flight path angle and heading angle, i.e., the target angular velocity, may be expressed as:

$\begin{matrix} {\overset{.}{\gamma} = \frac{a_{T}\cos\;\phi}{\upsilon_{T}}} & \lbrack 45\rbrack \\ {\overset{.}{\psi} = \frac{a_{T}\sin\;\phi}{\upsilon_{T}\cos\;\gamma}} & \; \end{matrix}$

Finally, there is a pseudo-measurement formed based on the orthogonality between the target velocity and acceleration vectors. z _(p1) =a _(tx)(u _(r) +v _(1x))+a _(Ty)(v _(r) +v _(Ty))+a _(Tz)(ω_(r) +v _(lz))+v _(p1)  [46] where v_(p1) is a zero mean, white Gaussian process with variance V_(p1) representing the uncertainty in the orthogonality between the target acceleration and velocity vectors. Also, there are ten additional pseudo-measurements formed because of the trigonometry on the target states. z _(p2) =a ₁ ²+a₂ ² −a ₉ ² +u _(p2) z _(p3) =a ₃ ² +a ₄ ² −a ₁₀ ² +u _(p3) z _(p4) =a ₅ ² +a ₆ ² −a ₁₁ ² +u _(p4) z _(p5) =a ₇ ² +a ₈ ² −a ₁₂ ² +u _(p5) z _(p6) =a ₁ ² +a ₃ ² −a ₁₃ ² +u _(p6) z _(p7) =a ₂ ² +a ₄ ² −a ₁₄ ² +u _(p7) z _(p8) =a ₅ ² +a ₇ ² −a ₁₅ ² +u _(p8) z _(p9) =a ₆ ² +a ₈ ² −a ₁₆ ² +u _(p9) z _(p10) =a ₉ ² +a ₁₀ ² −a ₁₃ ² −a ₁₄ ² +u _(p10) z _(p11) =a ₁₁ ² +a ₁₂ ² −a ₁₅ ² −a ₁₆ ² +u _(p11)  [47] where v_(p2), v_(p3) . . . v_(p11) are zero mean, white, Gaussian processes with variances V_(p2), V_(p3) . . . V_(p11), respectively, representing the uncertainty in the target states. When these pseudo-measurements are implemented, they are also set to zero.

3.2 Target Acceleration Estimator

In this section, the TAE is developed by using the extended Kalman filter to estimate the engagement states based on equation 33.1 to equation 33.6 and equation 44. By assuming the interceptor has measurements of its own acceleration (i.e., ā_(1x), ā_(1y) and ā_(1z).), equation 33.7 to equation 33.9 are not needed in the filter design. The states x, inputs u, state-dependent process noise ω, non-state-dependent process noise {tilde over (w)}, measurements z and sensor noise u may be represented with the following set of equations: x=[x _(r) y _(r) z _(r) u _(r) u _(r) ω_(r) a ₁ a ₂ . . . a ₁₆ φ]^(T) u=[ā _(1x) ā_(1y) ā_(1z)]^(T) ω=[{dot over (ω)}_(γ) {dot over (ω)}_(ψ)]^(T) {tilde over (ω)}={dot over (ω)}_(φ)  [48] z=[z_(el) z_(az) z_(R) z_({dot over (R)}) z_(p1) z_(p2) . . . z_(p11)]^(T) v=[v_(el) v_(az) v_(R) v_({dot over (R)}) v_(p1) v_(p2) . . . v_(p11)]^(T)

Accordingly, the system dynamics and measurement equation for the filter design can be written {dot over (x)}=f(x)+Bu+G(x)ω+{tilde over (G)}{tilde over (ω)}  [49] as z=h(x)+v where f(x) may be represented by equations 33.1 to equation 33.6 and equation 44, and h(x) may be represented by equation 34, equation 46 and equation 47.

Since the system dynamics and measurement equation are nonlinear, the extended Kalman filter, which includes a time propagation step and a measurement update step, may be used to estimate the states. The a priori state estimate may be represented as x(t|t_(i−1)) and the a priori error covariance may be denoted as P(t|t_(i−1)) where t_(i−1≦t)≦t_(i). The a posteriori state estimate may be denoted as {circumflex over (x)}(t_(i)) and a posteriori error covariance may be denoted as P(t_(i)). The state covariance E[x(t)x(t)^(T)] may be represented as X(t) because of the state-dependent process noise. For the time propagation, x(t|t_(i−1)) and P(t|t_(i−1)) are propagated from the current time t_(i−1) to the next sample time t_(i) by integrating {dot over (x)}( t|t _(i−1))=f[ x (t|t _(i−1))]+Bu(t) {dot over (P)}( t|t _(i−1))=A(t) P (t|t _(i−1))A(t)^(T) +E{G[x(t)]ΩG[x(t)]^(T) }+{tilde over (G)}{tilde over (Ω)}{tilde over (G)} ^(T)  [50] with the initial condition x(t_(i−1)|t_(i−1))={circumflex over (x)}(t_(i−1)) and P(t_(i−1)|t_(i−1))=P(t_(i−1)) where

$\begin{matrix} {{A(t)} = {\frac{\partial f}{\partial x}❘_{x = {\overset{\_}{x}{({t❘t_{i - 1}})}}}}} & \lbrack 51\rbrack \end{matrix}$

The term due to the state-dependent process noise E{G[x(t)]ΩG[x(t)]^(T)} may be determined by integrating {dot over (X)}(t)=A(t)X(t)+X(t)A(t)^(T) +E{G[x(t)]ΩG[x(t)]^(T) }+{tilde over (G)}{tilde over (Ω)}{tilde over (G)} ^(T)  [52] where

$\begin{matrix} {\Omega = \begin{bmatrix} \Omega_{\gamma} & 0 \\ 0 & \Omega_{\psi} \end{bmatrix}} & \lbrack 53\rbrack \\ {\overset{\sim}{\Omega} = \Omega_{\phi}} & \; \end{matrix}$

For the measurement update, {circumflex over (x)}(t_(i)) and P(t_(i)) may be updated by applying the set of equations: {circumflex over (x)}(t _(i))= x(t _(i) |t _(i−1))+K(t _(i)){z(t _(i))−h| x (t _(i) |t _(i−1))]} P(t _(i))=[I−K(t _(i))H(_(i))] P (t _(i) |t _(i−1))[I−K(t _(i))H(t _(i))]^(T) +K(t _(i))VK(t _(i))^(T)  [54] where the filter gain may be updated according to the following representation: K(t _(i))= P (t _(i) |t _(i−1))H(t _(i))^(t) [H(t _(i)) P (t _(i) |t _(i−1))H(t _(i))^(T) +V] ⁻¹  [55] and the linearized measurement matrix may be represented and updated according to the following

$\begin{matrix} {{H\left( t_{i} \right)} = {\frac{\partial h}{\partial x}❘_{x = {\overset{\_}{x}{({t❘t_{i - 1}})}}}}} & \lbrack 56\rbrack \end{matrix}$ and V is a diagonal matrix with V_(el), V_(az), V_(R), V_(R) and V_(p1), V_(p2) . . . V_(p11) on the diagonal line.

3.3 Linear Quadratic Guidance Law

In this section, the LQGL structure is explained by a solution to a deterministic linear quadratic optimization problem. For the derivation of the guidance law, it may be assumed that the target angular velocity {dot over (γ)} and {dot over (ψ)} are constant. The optimization problem for determining the guidance law may be expressed as:

$\begin{matrix} {{\min\limits_{a_{Ix},a_{Iy},a_{Iz}}{\frac{1}{2}\left\lbrack {{x_{r}\left( t_{f} \right)}^{2} + {y_{r}\left( t_{f} \right)}^{2} + {z_{r}\left( t_{f} \right)}^{2}} \right\rbrack}} + {\frac{c}{2}{\int_{t}^{t_{f}}{\left\lbrack {{a_{Ix}(\tau)}^{2} + {a_{Iy}(\tau)}^{2} + {a_{Iz}(\tau)}^{2}} \right\rbrack\ {\mathbb{d}\tau}}}}} & \lbrack 57\rbrack \end{matrix}$ subject to equation 33 and equation 44.1 to equation 44.3 without the noise {dot over (ω)}_(γ) and {dot over (ω)}_(γ) where c is a design parameter, t is the current time and t_(f) is the final time. The optimal solutions for the intercept acceleration command a*_(1x), a*_(1y) and a*_(1z) are

$\begin{matrix} {{a_{Ix}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 58.1\rbrack \\ \left\{ {{x_{r}(t)} + {T_{go}{u_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Ix}(t)}} - {{\left\lbrack {0\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu}\text{-1}\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\gamma\psi\phi}^{(2)}(t)} + {T_{go}{\Phi_{\gamma\psi\phi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma\psi\phi}\left( T_{go} \right)}\left\lbrack {{a_{1}(t)}\mspace{14mu}{a_{2}(t)}\mspace{14mu}{a_{3}(t)}\mspace{14mu}{a_{4}(t)}\mspace{14mu}{a_{5}(t)}\mspace{14mu}{a_{6}(t)}\mspace{14mu}{a_{7}(t)}\mspace{14mu}{a_{8}(t)}} \right\rbrack}^{T}} -} \right. & \; \\ \left. {{\left\lbrack {1\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\psi\phi}^{(2)}(t)} + {T_{go}{\Phi_{\psi\phi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\psi\phi}\left( T_{go} \right)}\left\lbrack {{a_{13}(t)}\mspace{14mu}{a_{14}(t)}\mspace{14mu}{a_{15}(t)}\mspace{14mu}{a_{16}(t)}} \right\rbrack}^{T}} \right\} & \; \\ {{a_{Iy}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 58.2\rbrack \\ \left\{ {{y_{r}(t)} + {T_{go}{\upsilon_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Iy}(t)}} - {{\left\lbrack {0\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 1} \right\rbrack\left\lbrack {{\Phi_{\gamma\phi}^{(2)}(t)} + {T_{go}{\Phi_{\gamma\phi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma\phi}\left( T_{go} \right)}\left\lbrack {{a_{9}(t)}\mspace{14mu}{a_{10}(t)}\mspace{14mu}{a_{11}(t)}\mspace{14mu}{a_{12}(t)}} \right\rbrack}^{T}}} \right\} & \; \\ {{a_{Iz}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 58.3\rbrack \\ \left\{ {{z_{r}(t)} + {T_{go}{\omega_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Iz}(t)}} - {{\left\lbrack {0\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0\mspace{20mu}\text{-1}\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\gamma\psi\phi}^{(2)}(t)} + {T_{go}{\Phi_{\gamma\psi\phi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma\psi\phi}\left( T_{go} \right)}\left\lbrack {{a_{1}(t)}\mspace{14mu}{a_{2}(t)}\mspace{14mu}{a_{3}(t)}\mspace{14mu}{a_{4}(t)}\mspace{14mu}{a_{5}(t)}\mspace{14mu}{a_{6}(t)}\mspace{14mu}{a_{7}(t)}\mspace{14mu}{a_{8}(t)}} \right\rbrack}^{T}} - {{\left\lbrack {0\mspace{14mu}\text{-1}\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\psi\phi}^{(2)}(t)} + {T_{go}{\Phi_{\psi\phi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\psi\phi}\left( T_{go} \right)}\left\lbrack {{a_{13}(t)}\mspace{14mu}{a_{14}(t)}\mspace{14mu}{a_{15}(t)}\mspace{14mu}{a_{16}(t)}} \right\rbrack}^{T}}} \right\} & \; \end{matrix}$ Φ_(γvφ), Φ_(γφ) and Φ_(ψφ) are defined as

$\begin{matrix} {{{{{{\Phi_{\gamma\psi\phi}(\tau)} =}\quad}\quad}\quad}{\quad{\frac{1}{4}{{\mathbb{e}}^{- \frac{{({\Omega_{\gamma} + \Omega_{\psi}})}\tau}{2}}\begin{bmatrix} {\Phi_{\gamma\psi\phi 1}(\tau)} & {\Phi_{\gamma\psi\phi 2}(\tau)} & {\Phi_{\gamma\psi\phi 3}(\tau)} & {\Phi_{\gamma\psi\phi 4}(\tau)} & \ldots & {\Phi_{\gamma\psi\phi 5}(\tau)} & {\Phi_{\gamma\psi\phi 6}(\tau)} & {\Phi_{\gamma\psi\phi 7}(\tau)} & {\Phi_{\gamma\psi\phi 8}(\tau)} \\ {- {\Phi_{\gamma\psi\phi 2}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} & {- {\Phi_{\gamma\psi\phi 4}(\tau)}} & {\Phi_{\gamma\psi\phi 3}(\tau)} & \ldots & {- {\Phi_{\gamma\psi\phi 6}(\tau)}} & {\Phi_{\gamma\psi\phi 5}(\tau)} & {- {\Phi_{\gamma\psi\phi 8}(\tau)}} & {\Phi_{\gamma\psi\phi 7}(\tau)} \\ {- {\Phi_{\gamma\psi\phi 3}(\tau)}} & {- {\Phi_{\gamma\psi\phi 4}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} & {\Phi_{\gamma\psi\phi 2}(\tau)} & \ldots & {- {\Phi_{\gamma\psi\phi 7}(\tau)}} & {- {\Phi_{\gamma\psi\phi 8}(\tau)}} & {\Phi_{\gamma\psi\phi 5}(\tau)} & {\Phi_{\gamma\psi\phi 6}(\tau)} \\ {\Phi_{\gamma\psi\phi 4}(\tau)} & {- {\Phi_{\gamma\psi\phi 3}(\tau)}} & {- {\Phi_{\gamma\psi\phi 2}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} & \ldots & {\Phi_{\gamma\psi\phi 8}(\tau)} & {- {\Phi_{\gamma\psi\phi 7}(\tau)}} & {- {\Phi_{\gamma\psi\phi 6}(\tau)}} & {\Phi_{\gamma\psi\phi 5}(\tau)} \\ {- {\Phi_{\gamma\psi\phi 5}(\tau)}} & {- {\Phi_{\gamma\psi\phi 6}(\tau)}} & {- {\Phi_{\gamma\psi\phi 7}(\tau)}} & {- {\Phi_{\gamma\psi\phi 8}(\tau)}} & \ldots & {\Phi_{\gamma\psi\phi 1}(\tau)} & {\Phi_{\gamma\psi\phi 2}(\tau)} & {\Phi_{\gamma\psi\phi 3}(\tau)} & {\Phi_{\gamma\psi\phi 4}(\tau)} \\ {\Phi_{\gamma\psi\phi 6}(\tau)} & {- {\Phi_{\gamma\psi\phi 5}(\tau)}} & {\Phi_{\gamma\psi\phi 8}(\tau)} & {- {\Phi_{\gamma\psi\phi 7}(\tau)}} & \ldots & {- {\Phi_{\gamma\psi\phi 2}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} & {- {\Phi_{\gamma\psi\phi 4}(\tau)}} & {\Phi_{\gamma\psi\phi 3}(\tau)} \\ {\Phi_{\gamma\psi\phi 7}(\tau)} & {\Phi_{\gamma\psi\phi 8}(\tau)} & {- {\Phi_{\gamma\psi\phi 5}(\tau)}} & {- {\Phi_{\gamma\psi\phi 6}(\tau)}} & \ldots & {- {\Phi_{\gamma\psi\phi 3}(\tau)}} & {- {\Phi_{\gamma\psi\phi 4}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} & {\Phi_{\gamma\psi\phi 2}(\tau)} \\ {- {\Phi_{\gamma\psi\phi 8}(\tau)}} & {\Phi_{\gamma\psi\phi 7}(\tau)} & {\Phi_{\gamma\psi\phi 6}(\tau)} & {- {\Phi_{\gamma\psi\phi 5}(\tau)}} & \ldots & {\Phi_{\gamma\psi\phi 4}(\tau)} & {- {\Phi_{\gamma\psi\phi 3}(\tau)}} & {- {\Phi_{\gamma\psi\phi 2}(\tau)}} & {\Phi_{\gamma\psi\phi 1}(\tau)} \end{bmatrix}}}}} & \lbrack 59\rbrack \end{matrix}$ where

$\begin{matrix} {{\Phi_{\gamma\psi\phi 1}(\tau)} = {{{\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {{\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {\ldots\mspace{14mu}{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \lbrack 60\rbrack \\ {{\Phi_{\gamma\psi\phi 2}(\tau)} = {{{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} - {{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} - {\ldots\mspace{14mu}{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 3}(\tau)} = {{{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {\ldots\mspace{14mu}{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 4}(\tau)} = {{{- {\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}}\tau} - {{\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {\ldots\mspace{14mu}{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 5}(\tau)} = {{{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} - {{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} - {\ldots\mspace{14mu}{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 6}(\tau)} = {{{- {\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}}\tau} + {{\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} - {\ldots\mspace{14mu}{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 7}(\tau)} = {{{- {\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}}\tau} + {{\cos\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} - {{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} + {\ldots\mspace{14mu}{\cos\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {{\Phi_{\gamma\psi\phi 8}(\tau)} = {{{- {\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}} \right)}}\tau} + {{\sin\left( {\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau} + {{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}} \right)}\tau} - {\ldots\mspace{14mu}{\sin\left( {\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}} \right)}\tau}}} & \; \\ {and} & \; \\ {{\Phi_{\gamma\phi}(\tau)} = {\frac{1}{2}{{\mathbb{e}}^{- \frac{\Omega_{\gamma}\tau}{2}}\begin{bmatrix} {\Phi_{\gamma\phi 1}(\tau)} & {\Phi_{\gamma\phi 2}(\tau)} & {\Phi_{\gamma\phi 3}(\tau)} & {\Phi_{\gamma\phi 4}(\tau)} \\ {- {\Phi_{\gamma\phi 2}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} & {- {\Phi_{\gamma\phi 4}(\tau)}} & {\Phi_{\gamma\phi 3}(\tau)} \\ {- {\Phi_{\gamma\phi 3}(\tau)}} & {- {\Phi_{\gamma\phi 4}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} & {\Phi_{\gamma\phi 2}(\tau)} \\ {\Phi_{\gamma\phi 4}(\tau)} & {- {\Phi_{\gamma\phi 3}(\tau)}} & {- {\Phi_{\gamma\phi 2}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} \end{bmatrix}}}} & \lbrack 61\rbrack \end{matrix}$ where φ_(γφ1)(τ)=cos ({dot over (γ)}+{dot over (φ)})τ+cos ({dot over (γ)}−{dot over (φ)})τ φ_(γφ2)(τ)=sin ({dot over (γ)}+{dot over (φ)})τ+sin ({dot over (γ)}−{dot over (φ)})τ φ_(γφ3)(τ)=sin ({dot over (γ)}+{dot over (φ)})τ−sin ({dot over (γ)}−{dot over (φ)})τ φ_(γφ4)(τ)=−cos ({dot over (γ)}+{dot over (φ)})τ+cos ({dot over (γ)}−{dot over (φ)})τ  [62] and

$\begin{matrix} {{\Phi_{\psi\phi}(\tau)} = {\frac{1}{2}{{\mathbb{e}}^{- \frac{\Omega_{\psi}\tau}{2}}\begin{bmatrix} {\Phi_{\psi\phi 1}(\tau)} & {\Phi_{\psi\phi 2}(\tau)} & {\Phi_{\psi\phi 3}(\tau)} & {\Phi_{\psi\phi 4}(\tau)} \\ {- {\Phi_{\psi\phi 2}(\tau)}} & {\Phi_{\psi\phi 1}(\tau)} & {- {\Phi_{\psi\phi 4}(\tau)}} & {\Phi_{\psi\phi 3}(\tau)} \\ {- {\Phi_{\psi\phi 3}(\tau)}} & {- {\Phi_{\psi\phi 4}(\tau)}} & {\Phi_{\psi\phi 1}(\tau)} & {\Phi_{\psi\phi 2}(\tau)} \\ {\Phi_{\psi\phi 4}(\tau)} & {- {\Phi_{\psi\phi 3}(\tau)}} & {- {\Phi_{\psi\phi 2}(\tau)}} & {\Phi_{\psi\phi 1}(\tau)} \end{bmatrix}}}} & \lbrack 63\rbrack \end{matrix}$ where φ_(ψφ1)(τ)=cos ({dot over (ψ)}+{dot over (φ)})τ+cos ({dot over (ψ)}−{dot over (φ)})τ φ_(ψφ2)(τ)=sin ({dot over (ψ)}+{dot over (φ)})τ+sin ({dot over (ψ)}−{dot over (φ)})τ  [64] φ_(ψφ3)(τ)=sin ({dot over (ψ)}+{dot over (φ)})τ−sin ({dot over (ψ)}−{dot over (φ)})τ φ_(ψφ4)(τ)=−cos ({dot over (ψ)}+{dot over (φ)})τ+cos ({dot over (ψ)}−{dot over (φ)})τ φ_(γψφ) ⁽¹⁾,φ_(γφ) ⁽¹⁾ and φ_(ψφ) ⁽¹⁾ are defined as

$\begin{matrix} {{\Phi_{\gamma\psi\phi}^{(1)}(t)} = {\frac{1}{4}\begin{bmatrix} {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 3}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 4}^{(1)}(t)} & \ldots & {\Phi_{\gamma\psi\phi 5}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 7}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 8}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 3}^{(1)}(t)} & \ldots & {- {\Phi_{\gamma\psi\phi 6}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 8}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 7}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(1)}(t)} & \ldots & {- {\Phi_{\gamma\psi\phi 7}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 8}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(1)}(t)} \\ {\Phi_{\gamma\psi\phi 4}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & \ldots & {\Phi_{\gamma\psi\phi 8}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 7}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi\phi 5}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 7}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 8}^{(1)}(t)}} & \ldots & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 3}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 4}^{(1)}(t)} \\ {\Phi_{\gamma\psi\phi 6}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 8}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 7}^{(1)}(t)}} & \ldots & {- {\Phi_{\gamma\psi\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 3}^{(1)}(t)} \\ {\Phi_{\gamma\psi\phi 7}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 8}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(1)}(t)}} & \ldots & {- {\Phi_{\gamma\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi\phi 8}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 7}^{(1)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(1)}(t)}} & \ldots & {\Phi_{\gamma\psi\phi 4}^{(1)}(t)} & {- {\Phi_{\gamma\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(1)}(t)} \end{bmatrix}}} & \lbrack 65\rbrack \\ {{\Phi_{\gamma\psi\phi 1}^{(1)}(t)} = {{\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \lbrack 66\rbrack \\ {{\Phi_{\gamma\psi\phi 2}^{(1)}(t)} = {{\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 3}^{(1)}(t)} = {{\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 4}^{(1)}(t)} = {{- {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} - {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 5}^{(1)}(t)} = {{\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 6}^{(1)}(t)} = {{- {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 7}^{(1)}(t)} = {{- {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} - {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 8}^{(1)}(t)} = {{- {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ \begin{matrix} {{\Phi_{\gamma\phi}^{(1)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\gamma\phi 1}^{(1)}(t)} & {\Phi_{\gamma\phi 2}^{(1)}(t)} & {\Phi_{\gamma\phi 3}^{(1)}(t)} & {\Phi_{\gamma\phi 4}^{(1)}(t)} \\ {- {\Phi_{\gamma\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\phi 1}^{(1)}(t)} & {- {\Phi_{\gamma\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\phi 3}^{(1)}(t)} \\ {- {\Phi_{\gamma\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\phi 4}^{(1)}(t)}} & {\Phi_{\gamma\phi 1}^{(1)}(t)} & {\Phi_{\gamma\phi 2}^{(1)}(t)} \\ {\Phi_{\gamma\phi 4}^{(1)}(t)} & {- {\Phi_{\gamma\phi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\phi 2}^{(1)}(t)}} & {\Phi_{\gamma\phi 1}^{(1)}(t)} \end{bmatrix}}} & \lbrack 67\rbrack \end{matrix} & \; \end{matrix}$ φ_(γφ1) ⁽¹⁾(t)=φ_(cos) ⁽¹⁾(Ω_(γ),{dot over (γ)}+{dot over (φ)})+φ_(cos) ⁽¹⁾(Ω_(γ),{dot over (γ)}−{dot over (φ)}) φ_(γφ2) ⁽¹⁾(t)=φ_(sin) ⁽¹⁾(Ω_(γ),{dot over (γ)}+{dot over (φ)})+φ_(sin) ⁽¹⁾(Ω_(γ), {dot over (γ)}−{dot over (φ)}) φ_(γφ3) ⁽¹⁾(t)=φ_(sin) ⁽¹⁾(Ω_(γ), {dot over (γ)}+{dot over (φ)})+φ_(sin) ⁽¹⁾(Ω_(γ), {dot over (γ)}−{dot over (φ)}) φ_(γφ4) ⁽¹⁾(t)=−φ_(cos) ⁽¹⁾(Ω_(γ),{dot over (γ)}+{dot over (φ)})+φ_(cos) ⁽¹⁾(Ω_(γ), {dot over (γ)}−{dot over (φ)})   [68]

$\begin{matrix} {{\Phi_{\psi\phi}^{(1)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\psi\phi 1}^{(1)}(t)} & {\Phi_{\psi\phi 2}^{(1)}(t)} & {\Phi_{\psi\phi 3}^{(1)}(t)} & {\Phi_{\psi\phi 4}^{(1)}(t)} \\ {- {\Phi_{\psi\phi 2}^{(1)}(t)}} & {\Phi_{\psi\phi 1}^{(1)}(t)} & {- {\Phi_{\psi\phi 4}^{(1)}(t)}} & {\Phi_{\psi\phi 3}^{(1)}(t)} \\ {- {\Phi_{\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\psi\phi 4}^{(1)}(t)}} & {\Phi_{\psi\phi 1}^{(1)}(t)} & {\Phi_{\psi\phi 2}^{(1)}(t)} \\ {\Phi_{\psi\phi 4}^{(1)}(t)} & {- {\Phi_{\psi\phi 3}^{(1)}(t)}} & {- {\Phi_{\psi\phi 2}^{(1)}(t)}} & {\Phi_{\psi\phi 1}^{(1)}(t)} \end{bmatrix}}} & \lbrack 69\rbrack \end{matrix}$ φ_(ψφ1) ⁽¹⁾(t)=φ_(cos) ⁽¹⁾(Ω_(ψ),{dot over (ψ)}+{dot over (φ)})+φ_(cos) ⁽¹⁾(Ω_(ψ),{dot over (ψ)}−{dot over (φ)}) φ_(ψφ2) ⁽¹⁾(t)=φ_(sin) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}+{dot over (φ)})+φ_(sin) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)}) φ_(ψφ3) ⁽¹⁾(t)=φ_(sin) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}+{dot over (φ)})−φ_(sin) ⁽¹⁾(Ω_(sin) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)}) φ_(ψφ4) ⁽¹⁾(t)=−φ_(cos) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}+{dot over (φ)})+φ_(cos) ⁽¹⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)})  [70] where

$\begin{matrix} {{{\Phi_{\sin}^{(1)}\left( {\Omega,\omega} \right)} = {\frac{1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\omega\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} + {\frac{\Omega}{2}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}{\Phi_{\cos}^{(1)}\left( {\Omega,\omega} \right)} = {\frac{1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\frac{\Omega}{2}\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} - {{\omega\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}} & \lbrack 71\rbrack \end{matrix}$ φ_(γψφ) ⁽²⁾, φ_(γφ) ⁽²⁾ and φ_(ψφ) ⁽²⁾ may be defined as

$\begin{matrix} {{\Phi_{\gamma\psi\phi}^{(2)}(t)} = {\frac{1}{4}\begin{bmatrix} {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 3}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 4}^{(2)}(t)} & \ldots & {\Phi_{\gamma\psi\phi 5}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 7}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 8}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 3}^{(2)}(t)} & \ldots & {- {\Phi_{\gamma\psi\phi 6}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 8}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 7}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(2)}(t)} & \ldots & {- {\Phi_{\gamma\psi\phi 7}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 8}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(2)}(t)} \\ {\Phi_{\gamma\psi\phi 4}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & \ldots & {\Phi_{\gamma\psi\phi 8}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 7}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 5}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi\phi 5}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 7}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 8}^{(2)}(t)}} & \ldots & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 3}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 4}^{(2)}(t)} \\ {\Phi_{\gamma\psi\phi 6}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 8}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 7}^{(2)}(t)}} & \ldots & {- {\Phi_{\gamma\psi\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 3}^{(2)}(t)} \\ {\Phi_{\gamma\psi\phi 7}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 8}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 6}^{(2)}(t)}} & \ldots & {- {\Phi_{\gamma\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 2}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi\phi 8}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 7}^{(2)}(t)} & {\Phi_{\gamma\psi\phi 6}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 5}^{(2)}(t)}} & \ldots & {\Phi_{\gamma\psi\phi 4}^{(2)}(t)} & {- {\Phi_{\gamma\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi\phi 1}^{(2)}(t)} \end{bmatrix}}} & \lbrack 72\rbrack \\ {{\Phi_{\gamma\psi\phi 1}^{(2)}(t)} = {{\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \lbrack 73\rbrack \\ {{\Phi_{\gamma\psi\phi 2}^{(2)}(t)} = {{\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 3}^{(2)}(t)} = {{\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 4}^{(2)}(t)} = {{- {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} - {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 5}^{(2)}(t)} = {{\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 6}^{(2)}(t)} = {{- {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 7}^{(2)}(t)} = {{- {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} - {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ {{\Phi_{\gamma\psi\phi 8}^{(2)}(t)} = {{- {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} + \overset{.}{\phi}}} \right)}} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi} - \overset{.}{\phi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} + \overset{.}{\phi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi} - \overset{.}{\phi}}} \right)}}} & \; \\ \begin{matrix} {{\Phi_{\gamma\phi}^{(2)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\gamma\phi 1}^{(2)}(t)} & {\Phi_{\gamma\phi 2}^{(2)}(t)} & {\Phi_{\gamma\phi 3}^{(2)}(t)} & {\Phi_{\gamma\phi 4}^{(2)}(t)} \\ {- {\Phi_{\gamma\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\phi 1}^{(2)}(t)} & {- {\Phi_{\gamma\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\phi 3}^{(2)}(t)} \\ {- {\Phi_{\gamma\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\phi 4}^{(2)}(t)}} & {\Phi_{\gamma\phi 1}^{(2)}(t)} & {\Phi_{\gamma\phi 2}^{(2)}(t)} \\ {\Phi_{\gamma\phi 4}^{(2)}(t)} & {- {\Phi_{\gamma\phi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\phi 2}^{(2)}(t)}} & {\Phi_{\gamma\phi 1}^{(2)}(t)} \end{bmatrix}}} & \lbrack 74\rbrack \end{matrix} & \; \end{matrix}$ φ_(γφ1) ⁽²⁾(t)=φ_(cos) ⁽²⁾(Ω_(γ), {dot over (γ)}+{dot over (φ)})+φ_(cos) ⁽²⁾(Ω_(γ),{dot over (γ)}−{dot over (φ)}) φ_(γφ2) ⁽²⁾(t)=φ_(sin) ⁽²⁾(Ω_(γ), {dot over (γ)}+{dot over (φ)})+φ_(sin) ⁽²⁾(Ω_(γ),{dot over (γ)}−{dot over (φ)}) φ_(γφ3) ⁽²⁾(t)=φ_(sin) ⁽²⁾(Ω_(γ), {dot over (γ)}+{dot over (φ)})−φ_(sin) ⁽²⁾(Ω_(γ),{dot over (γ)}−{dot over (φ)}) φ_(γφ4) ⁽²⁾(t)=−φ_(cos) ⁽²⁾(Ω_(γ), {dot over (γ)}+{dot over (φ)})+_(100 cos) ⁽²⁾(Ω_(γ), {dot over (γ)}−{dot over (φ)})  [75]

$\begin{matrix} {{\Phi_{\psi\phi}^{(2)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\psi\phi 1}^{(2)}(t)} & {\Phi_{\psi\phi 2}^{(2)}(t)} & {\Phi_{\psi\phi 3}^{(2)}(t)} & {\Phi_{\psi\phi 4}^{(2)}(t)} \\ {- {\Phi_{\psi\phi 2}^{(2)}(t)}} & {\Phi_{\psi\phi 1}^{(2)}(t)} & {- {\Phi_{\psi\phi 4}^{(2)}(t)}} & {\Phi_{\psi\phi 3}^{(2)}(t)} \\ {- {\Phi_{\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\psi\phi 4}^{(2)}(t)}} & {\Phi_{\psi\phi 1}^{(2)}(t)} & {\Phi_{\psi\phi 2}^{(2)}(t)} \\ {\Phi_{\psi\phi 4}^{(2)}(t)} & {- {\Phi_{\psi\phi 3}^{(2)}(t)}} & {- {\Phi_{\psi\phi 2}^{(2)}(t)}} & {\Phi_{\psi\phi 1}^{(2)}(t)} \end{bmatrix}}} & \lbrack 76\rbrack \end{matrix}$ φ_(ψφ1) ⁽²⁾(t)=φ_(cos) ⁽²⁾(Ω_(ψ),{dot over (ψ)}+{dot over (φ)})+φ_(cos) ⁽²⁾(Ω_(ψ),{dot over (ψ)}−{dot over (φ)}) φ_(ψφ2) ⁽²⁾(t)=φ_(sin) ⁽²⁾(Ω_(ψ),{dot over (ψ)}+{dot over (φ)})+φ_(sin) ⁽²⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)}) φ_(ψφ3)(t)=φ_(sin) ⁽²⁾(Ω_(ψ), {dot over (ψ)}+{dot over (φ)})−φ_(sin) ⁽²⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)})  [77] φ_(ψφ4) ⁽²⁾(t)=−φ_(cos) ⁽²⁾(Ω_(ψ), {dot over (ψ)}+{dot over (φ)})+φ_(cos) ⁽²⁾(Ω_(ψ), {dot over (ψ)}−{dot over (φ)}) where

$\begin{matrix} {{{\Phi_{\sin}^{(2)}\left( {\Omega,\omega} \right)} = {\frac{- 1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\omega\; T_{go}} + {\frac{\omega\Omega}{\frac{\Omega^{2}}{4} + \omega^{2}}\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} + {\frac{\frac{\Omega^{2}}{4} - \omega^{2}}{\frac{\Omega^{2}}{4} + \omega^{2}}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}{{\Phi_{\cos}^{(2)}\left( {\Omega,\omega} \right)} = {\frac{- 1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {\frac{\Omega\; T_{go}}{2} + {\frac{\frac{\Omega^{2}}{4} - \omega^{2}}{\frac{\Omega^{2}}{4} + \omega^{2}}\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} - {\frac{\omega\Omega}{\frac{\Omega^{2}}{4} + \omega^{2}}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}} & \lbrack 78\rbrack \end{matrix}$

These solutions are functions of x_(r), y_(r), z_(r), u_(r), v_(r), ω_(r), ā_(1x), ā_(1y), ā_(1z), a₁,a₂ . . . a₁₆, {dot over (φ)}, {dot over (γ)}, {dot over (ψ)} and the time-to-go T_(go)=t_(f)−t. When implementing the guidance law, x_(r), y_(r), z_(r), u_(r), v_(r), ω_(r), a₁,a₂ . . . a₁₆ and {dot over (φ)} are estimated as {circumflex over (x)}. The intercept acceleration, ā_(1x), ā_(1y) and ā_(1z), are measured. The target angular velocity can be estimated by

$\begin{matrix} {{\overset{.}{\gamma} = \frac{{\overset{\rightarrow}{a}}_{T}{\cdot {\overset{\rightarrow}{a}}_{1}}}{\upsilon_{T}}}{\overset{.}{\psi} = \frac{{\overset{\rightarrow}{a}}_{T}{\cdot {\overset{\rightarrow}{a}}_{2}}}{\upsilon_{T}\cos\;\gamma}}} & \lbrack 79\rbrack \end{matrix}$ where v _(i)=√{square root over ((u _(r) +v _(1x))²+(v _(r) +v _(1y))²+(ω_(r) +v _(1z))²)}{square root over ((u _(r) +v _(1x))²+(v _(r) +v _(1y))²+(ω_(r) +v _(1z))²)}{square root over ((u _(r) +v _(1x))²+(v _(r) +v _(1y))²+(ω_(r) +v _(1z))²)}  [80] and ā_(T), ā₁ and ā₂ can be obtained from equation 42, equation 36 and equation 37, respectively, with

$\begin{matrix} {{\gamma = {\sin^{- 1}\frac{\upsilon_{r} + \upsilon_{Iy}}{\sqrt{\left( {u_{r} + \upsilon_{Ix}} \right)^{2} + \left( {\upsilon_{r} + \upsilon_{Iy}} \right)^{2} + \left( {\omega_{r} + \upsilon_{Iz}} \right)^{2}}}}}{\psi = {\tan^{- 1}\frac{\omega_{r} + \upsilon_{Iz}}{u_{r} + \upsilon_{Ix}}}}} & \lbrack 81\rbrack \end{matrix}$

The time-to-go can be estimated by

$\begin{matrix} {T_{go} = {\frac{R}{\overset{.}{R}} = \frac{x_{r}^{2} + y_{r}^{2} + z_{r}^{2}}{{x_{r}u_{r}} + {y_{r}\upsilon_{r}} + {z_{r}\omega_{r}}}}} & \lbrack 82\rbrack \end{matrix}$

3.4 Special Case

In this section, a special case of the TAE/LQGL algorithm is developed when the bank angle rate of the target acceleration vector is zero, i.e., the bank angle is a constant.

3.4.1 Target Model

When the bank angle is a constant, the target states are defined as a ₁ =a _(T) cos φ sin (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₂ =a _(T) cos φ sin (γ+ω_(γ)) cos (ψ+ω_(ψ)) a ₃ =a _(T) cos φ cos (γ+ω_(γ)) sin (ψ+ω_(ψ)) a ₄ =a _(T) cos φ cos (γ+ω_(γ)) cos (ψ+ω_(ψ))  [83] a ₅ =a _(T) cos φ sin (γ+ω_(γ)) a ₆ =a _(T) cos φ cos (γ+ω_(γ)) a ₇ =a _(T) sin φ sin (ψ+ω_(ψ)) a ₈ =a _(T) sin φ cos (ψ+ω_(ψ))

Then, equation 40 becomes a =(−a ₂ +a ₇)ī+a₆ j+(−a ₁ −a ₈) k   [84] and therefore a _(Tx) =−a ₂ +a ₇ a _(Ty) =a ₆ a _(Tz) =−a ₁ −a ₈  [85]

By taking the Ito stochastic differential of equation 83, the target model is obtained as

$\begin{matrix} {\begin{bmatrix} {\overset{.}{a}}_{1} \\ {\overset{.}{a}}_{2} \\ {\overset{.}{a}}_{3} \\ {\overset{.}{a}}_{4} \end{bmatrix} = {{\begin{bmatrix} {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \overset{.}{\psi} & \overset{.}{\gamma} & 0 \\ {- \overset{.}{\psi}} & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & 0 & \overset{.}{\gamma} \\ {- \overset{.}{\gamma}} & 0 & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} & \overset{.}{\psi} \\ 0 & {- \overset{.}{\gamma}} & {- \overset{.}{\psi}} & {- \frac{\Omega_{\gamma} + \Omega_{\psi}}{2}} \end{bmatrix}\begin{bmatrix} a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \end{bmatrix}} + {\begin{bmatrix} a_{3} & a_{2} \\ a_{4} & {- a_{1}} \\ {- a_{1}} & a_{4} \\ {- a_{2}} & {- a_{3}} \end{bmatrix}\begin{bmatrix} {\overset{.}{\omega}}_{\gamma} \\ {\overset{.}{\omega}}_{\psi} \end{bmatrix}}}} & \lbrack 86.1\rbrack \\ {\begin{bmatrix} {\overset{.}{a}}_{5} \\ {\overset{.}{a}}_{6} \end{bmatrix} = {{\begin{bmatrix} {- \frac{\Omega_{\gamma}}{2}} & \overset{.}{\gamma} \\ {- \overset{.}{\gamma}} & {- \frac{\Omega_{\gamma}}{2}} \end{bmatrix}\begin{bmatrix} a_{5} \\ a_{6} \end{bmatrix}} + {\begin{bmatrix} a_{6} \\ {- a_{5}} \end{bmatrix}{\overset{.}{\omega}}_{\gamma}}}} & \lbrack 86.2\rbrack \\ {\begin{bmatrix} {\overset{.}{a}}_{7} \\ {\overset{.}{a}}_{8} \end{bmatrix} = {{\begin{bmatrix} {- \frac{\Omega_{\psi}}{2}} & \overset{.}{\psi} \\ {- \overset{.}{\psi}} & {- \frac{\Omega_{\psi}}{2}} \end{bmatrix}\begin{bmatrix} a_{7} \\ a_{8} \end{bmatrix}} + {\begin{bmatrix} a_{8} \\ {- a_{7}} \end{bmatrix}{\overset{.}{\omega}}_{\psi}}}} & \lbrack 86.3\rbrack \end{matrix}$ where {dot over (ω)}_(γ) and {dot over (ω)}_(ψ) are zero mean, white Gaussian process with power spectral densities Ω_(γ) and Ω_(ψ), respectively, Finally, there are only two pseudo-measurements in the exemplary structural development because of the trigonometry on fit target states. z _(p2) =a ₁ ² +a ₂ ² −a ₅ ² +v _(p2) z _(p3) =a ₃ ² +a ₄ ² −a ₆ ² +v _(p3)  [87] where v_(p2) and v_(p3) are zero mean, white Gaussian process with variance v_(p2) and v_(p3), respectively, representing the uncertainty in the target states

3.4.2 Target Acceleration Estimator

In this section, the TAE structure is explained by using the extended Kalman filter to estimate the engagement states based on equation 33.1 to equation 33.6 and equation 86. Let the states x, inputs u, state-dependent process noise ω, measurements z and sensor noise v be x=[x_(r) y_(r) z_(r) u_(r) v_(r) ω_(r) a₁ a₂ . . . a₈]^(T) u=[ā_(1x) ā_(1y) ā_(1z)]^(T) ω=[{dot over (ω)}_(γ) {dot over (ω)}_(γ)]^(T)  [88] z=[z_(al) z_(az) z_(R) z_({dot over (R)}) z_(p1) z_(p2) z_(p3)]^(T) v=[v_(cl) v_(az) v_(R) v_({dot over (R)}) v_(p1) v_(p2) v_(p3)]^(T)

Then, the system dynamics and measurement equation for the filter design can be written as {dot over (x)}=f(x)+Bu+G(x)ω z=h(x)+v  [89] where f(x) represented by equation 33.1 to equation 33.6 and equation 86, and h(x) is represented by equation 34, equation 46 and equation 87.

Since the system dynamics and measurement equation are nonlinear, the extended Kalman filter, which includes a time propagation step and a measurement update step, may be used to estimate the states. The a priori state estimate may be represented as x(t═t_(i−1)) and the a priori error covariance may be represented as P(t|t_(i−1)) where t_(i−1)≦t≦t_(i). The a posteriori state estimate may be denoted as {circumflex over (x)}(t_(i)) and the a posteriori error covariance may be denoted as P(t_(i)). The state covariance E[x(t)x(t)^(T)] may be denoted as X(t) because of the state-dependent process noise. For the time propagation, x(t|t_(i−1)) and P(t|t_(i−1)) may be propagated from the current time t_(i−1) to the next sample time t_(i) by integrating {dot over (x)}( t|t _(i−1))=f[ x( t|t _(i−1))]+Bu(t) {dot over (P)}( t _(i−1))=A(t) P (t|t _(i−1))+ P (t|t _(i−1))A(t)^(T) +E{G[x(t)]ΩG[x(t)]^(T)}  [90] with the initial condition x(t_(i−1)|t_(i−1))={circumflex over (x)}(t_(i−1)) and P(t_(i−1)|t_(i−1))=P(t_(i−1)) where

$\begin{matrix} {{A(t)} = {\frac{\partial f}{\partial x}❘_{x = {\overset{\_}{x}{({t❘t_{i - 1}})}}}}} & \lbrack 91\rbrack \end{matrix}$

The term due to the state-dependent process noise E{G[x(t)]ΩG[x(t)]⁷} may be determined by integrating {dot over (X)}( t)=A(t)X(t)+X(t)A(t)^(T) +E{G[x(t)]ΩG[x(t)]^(T)}  [92] where

$\begin{matrix} {\Omega = \begin{bmatrix} \Omega_{\gamma} & 0 \\ 0 & \Omega_{\psi} \end{bmatrix}} & \lbrack 93\rbrack \end{matrix}$

For the measurement update, {circumflex over (x)}(t_(i)) and P(t_(i)) may be updated by applying {circumflex over (x)}(t _(i))= x(t _(i) |t _(i−1))+K(t _(i)){z(t _(i))−h[ x (t _(i) |t _(i−1))]} P(t _(i))=[I−K(t _(i))H(t _(i))] P (t _(i) |t _(i−1))[I−K(t _(i))H(t _(i))]^(T) +K(t _(i))VK(t _(i))^(T)  [94] where the filter gain may be represented as K(t _(i))= P (t _(i) |t _(i−1))H(t _(i))^(T) [H(t _(i)) P (t _(i) |t _(i−1))H(t _(i))^(T) +V] ⁻¹  [95] and the linearized measurement matrix may be represented as

$\begin{matrix} {{H\left( t_{i} \right)} = {\frac{\partial h}{\partial x}❘_{x = {\overset{\_}{x}{({t❘t_{i - 1}})}}}}} & \lbrack 96\rbrack \end{matrix}$ and V is a diagonal matrix with V_(el), V_(az), V_(R), V_({dot over (R)}), V_(p1), V_(p2) and V_(p3) on the diagonal line.

3.4.3 Linear Quadratic Guidance Law

In this section, the LQGL structure is described by solving a deterministic linear quadratic optimization problem. For the derivation of the guidance law, i.e., the LQGL, it is assumed that the heading angular velocity, {dot over (ψ)}, is constant. The optimization problem for determining the guidance law may be represented as

$\begin{matrix} {{\min\limits_{a_{Ix},a_{Iy},a_{Iz}}{{\frac{1}{2}\left\lbrack \quad \right.}{x_{r}\left( t_{f} \right)}^{2}}} + {y_{r}\left( t_{f} \right)}^{2} + {z_{r}\left( t_{f} \right)}^{2} + {\frac{c}{2}{\int_{t}^{t_{f}}{\left\lbrack {{a_{Ix}(\tau)}^{2} + {a_{I_{y}}(\tau)}^{2} + {a_{Iz}(\tau)}^{2}} \right\rbrack\ {\mathbb{d}\tau}}}}} & \lbrack 97\rbrack \end{matrix}$ subject to equation 33 and equation 86 without the noise {dot over (ω)}_(γ) and {dot over (ω)}_(ψ) where c is a design parameter, t is the current time and t_(f) is the final time. The optimal solutions of the above exemplary objective function of equation 97 for the intercept acceleration commands a*_(1x), a*_(1y) and a*_(1z) are as follows.

$\begin{matrix} {{a_{Ix}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 98.1\rbrack \\ \left\{ {{x_{r}(t)} + {T_{go}{u_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Ix}(t)}} - {{\left\lbrack {0\mspace{14mu}\text{-1}\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\gamma\psi}^{(2)}(t)} + {T_{go}{\Phi_{\gamma\psi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma\psi}\left( T_{go} \right)}\left\lbrack {{a_{1}(t)}\mspace{14mu}{a_{2}(t)}\mspace{14mu}{a_{3}(t)}\mspace{14mu}{a_{4}(t)}} \right\rbrack}^{T}} -} \right. & \; \\ \left. {{\left\lbrack {1\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\psi}^{(2)}(t)} + {T_{go}{\Phi_{\psi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\psi}\left( T_{go} \right)}\left\lbrack {{a_{7}(t)}\mspace{14mu}{a_{8}(t)}} \right\rbrack}^{T}} \right\} & \; \\ {{a_{Iy}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 98.2\rbrack \\ \left\{ {{y_{r}(t)} + {T_{go}{\upsilon_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Iy}(t)}} - {{\left\lbrack {0\mspace{14mu} 1} \right\rbrack\left\lbrack {{\Phi_{\gamma}^{(2)}(t)} + {T_{go}{\Phi_{\gamma}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma}\left( T_{go} \right)}\left\lbrack {{a_{5}(t)}\mspace{14mu}{a_{6}(t)}} \right\rbrack}^{T}}} \right\} & \; \\ {{a_{Iz}^{*}(t)} = {{- \frac{2}{\alpha}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right)\left( {{2\; c} + \frac{2T_{go}^{3}}{3} - \frac{2T_{go}^{2}}{\alpha} + \frac{2T_{go}}{\alpha^{2}} + \frac{1}{\alpha^{3}} - \frac{4T_{go}{\mathbb{e}}^{{- a}\; T_{go}}}{\alpha^{2}} - \frac{{\mathbb{e}}^{{- 2}\; a\; T_{go}}}{\alpha^{3}}} \right)^{- 1}}} & \lbrack 98.3\rbrack \\ \left\{ {{z_{r}(t)} + {T_{go}{\omega_{r}(t)}} + {\frac{1}{\alpha^{2}}\left( {1 - {\alpha\; T_{go}} - {\mathbb{e}}^{{- \alpha}\; T_{go}}} \right){{\overset{\_}{\alpha}}_{Iz}(t)}} - {{\left\lbrack {\text{-1}\mspace{14mu} 0\mspace{14mu} 0\mspace{14mu} 0} \right\rbrack\left\lbrack {{\Phi_{\gamma\psi}^{(2)}(t)} + {T_{go}{\Phi_{\gamma\psi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\gamma\psi}\left( T_{go} \right)}\left\lbrack {{a_{1}(t)}\mspace{14mu}{a_{2}(t)}\mspace{14mu}{a_{3}(t)}\mspace{14mu}{a_{4}(t)}} \right\rbrack}^{T}} - {{\left\lbrack {0\mspace{14mu}\text{-1}} \right\rbrack\left\lbrack {{\Phi_{\psi}^{(2)}(t)} + {T_{go}{\Phi_{\psi}^{(1)}(t)}}} \right\rbrack}{{\Phi_{\psi}\left( T_{go} \right)}\left\lbrack {{a_{7}(t)}\mspace{14mu}{a_{8}(t)}} \right\rbrack}^{T}}} \right\} & \; \end{matrix}$ φ_(γψ), φ_(γ) and φ_(ψ) are defined as

$\begin{matrix} {{\Phi_{\gamma\psi}(\tau)} = {\frac{1}{2}{{\mathbb{e}}^{- \frac{{({\Omega_{\gamma} + \Omega_{\psi}})}\tau}{2}}\begin{bmatrix} {\Phi_{\gamma\phi 1}(\tau)} & {\Phi_{\gamma\phi 2}(\tau)} & {\Phi_{\gamma\phi 3}(\tau)} & {\Phi_{\gamma\phi 4}(\tau)} \\ {- {\Phi_{\gamma\phi 2}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} & {- {\Phi_{\gamma\phi 4}(\tau)}} & {\Phi_{\gamma\phi 3}(\tau)} \\ {- {\Phi_{\gamma\phi 3}(\tau)}} & {- {\Phi_{\gamma\phi 4}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} & {\Phi_{\gamma\phi 2}(\tau)} \\ {\Phi_{\gamma\phi 4}(\tau)} & {- {\Phi_{\gamma\phi 3}(\tau)}} & {- {\Phi_{\gamma\phi 2}(\tau)}} & {\Phi_{\gamma\phi 1}(\tau)} \end{bmatrix}}}} & \lbrack 99\rbrack \end{matrix}$ where φ_(γψ1)(τ)=cos ({dot over (γ)}+{dot over (ψ)})τ+cos ({dot over (γ)}−{dot over (ψ)})τ φ_(γψ2)(τ)=sin ({dot over (γ)}+{dot over (ψ)})τ+sin ({dot over (γ)}−{dot over (ψ)})τ  [100] φ_(γψ3)(τ)=sin ({dot over (γ)}+{dot over (ψ)})τ+sin ({dot over (γ)}−{dot over (ψ)})τ φ_(γψ4)(τ)=−cos ({dot over (γ)}+{dot over (ψ)})τ+cos ({dot over (γ)}−{dot over (ψ)})τ and

$\begin{matrix} {{{\Phi_{\gamma}(\tau)} = {{\mathbb{e}}^{- \frac{\Omega_{\gamma}\tau}{2}}\begin{bmatrix} {\cos\;\overset{.}{\gamma}\tau} & {\sin\;\overset{.}{\gamma}\tau} \\ {{- \sin}\;\overset{.}{\gamma}\tau} & {\cos\;\overset{.}{\gamma}\tau} \end{bmatrix}}}{{\Phi_{\psi}(\tau)} = {{\mathbb{e}}^{- \frac{\Omega_{\psi}\tau}{2}}\begin{bmatrix} {\cos\;\overset{.}{\psi}\tau} & {\sin\;\overset{.}{\psi}\tau} \\ {{- \sin}\;\overset{.}{\psi}\tau} & {\cos\;\overset{.}{\psi}\tau} \end{bmatrix}}}} & \lbrack 101\rbrack \end{matrix}$ φ_(γψ) ⁽¹⁾, φ_(γ) ⁽¹⁾ and φ_(ψ) ⁽¹⁾ are defined as

$\begin{matrix} {{\Phi_{\gamma\psi}^{(1)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\gamma\psi 1}^{(1)}(t)} & {\Phi_{\gamma\psi 2}^{(1)}(t)} & {\Phi_{\gamma\psi 3}^{(1)}(t)} & {\Phi_{\gamma\psi 4}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi 1}^{(1)}(t)} & {- {\Phi_{\gamma\psi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi 3}^{(1)}(t)} \\ {- {\Phi_{\gamma\psi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi 4}^{(1)}(t)}} & {\Phi_{\gamma\psi 1}^{(1)}(t)} & {\Phi_{\gamma\psi 2}^{(1)}(t)} \\ {\Phi_{\gamma\psi 4}^{(1)}(t)} & {- {\Phi_{\gamma\psi 3}^{(1)}(t)}} & {- {\Phi_{\gamma\psi 2}^{(1)}(t)}} & {\Phi_{\gamma\psi 1}^{(1)}(t)} \end{bmatrix}}} & \lbrack 102\rbrack \\ \begin{matrix} {{\Phi_{\gamma\psi 1}^{(1)}(t)} = {{\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 2}^{(1)}(t)} = {{\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} - {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 3}^{(1)}(t)} = {{\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} + {\Phi_{\sin}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 4}^{(1)}(t)} = {{- {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)}} + {\Phi_{\cos}^{(1)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \end{matrix} & \lbrack 103\rbrack \\ {{\Phi_{\gamma}^{(1)}(t)} = \begin{bmatrix} {\Phi_{\cos}^{(1)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} & {\Phi_{\sin}^{(1)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} \\ {- {\Phi_{\sin}^{(1)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)}} & {\Phi_{\cos}^{(1)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} \end{bmatrix}} & \lbrack 104\rbrack \\ {{\Phi_{\gamma}^{(1)}(t)} = \begin{bmatrix} {\Phi_{\cos}^{(1)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} & {\Phi_{\sin}^{(1)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} \\ {- {\Phi_{\sin}^{(1)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)}} & {\Phi_{\cos}^{(1)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} \end{bmatrix}} & \lbrack 105\rbrack \end{matrix}$ where

$\begin{matrix} \begin{matrix} {{{\Phi_{\sin}^{(1)}\left( {\Omega,\omega} \right)} = {\frac{1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\omega\;\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} + {\frac{\Omega}{2}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}{\Phi_{\cos}^{(1)}\left( {\Omega,\omega} \right)} = {\frac{1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\frac{\Omega}{2}\;\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} - {{\omega\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}} & \lbrack 106\rbrack \end{matrix} & \; \end{matrix}$ φ_(γψ) ⁽²⁾, φ_(γ) ⁽²⁾ and φ_(ψ) ⁽²⁾ are defined as

$\begin{matrix} {{\Phi_{\gamma\psi}^{(2)}(t)} = {\frac{1}{2}\begin{bmatrix} {\Phi_{\gamma\psi 1}^{(2)}(t)} & {\Phi_{\gamma\psi 2}^{(2)}(t)} & {\Phi_{\gamma\psi 3}^{(2)}(t)} & {\Phi_{\gamma\psi 4}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi 1}^{(2)}(t)} & {- {\Phi_{\gamma\psi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi 3}^{(2)}(t)} \\ {- {\Phi_{\gamma\psi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi 4}^{(2)}(t)}} & {\Phi_{\gamma\psi 1}^{(2)}(t)} & {\Phi_{\gamma\psi 2}^{(2)}(t)} \\ {\Phi_{\gamma\psi 4}^{(2)}(t)} & {- {\Phi_{\gamma\psi 3}^{(2)}(t)}} & {- {\Phi_{\gamma\psi 2}^{(2)}(t)}} & {\Phi_{\gamma\psi 1}^{(2)}(t)} \end{bmatrix}}} & \lbrack 107\rbrack \\ \begin{matrix} {{\Phi_{\gamma\psi 1}^{(2)}(t)} = {{\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 2}^{(2)}(t)} = {{\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} - {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 3}^{(2)}(t)} = {{\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)} + {\Phi_{\sin}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \\ {{\Phi_{\gamma\psi 4}^{(2)}(t)} = {{- {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} + \overset{.}{\psi}}} \right)}} + {\Phi_{\cos}^{(2)}\left( {{\Omega_{\gamma} + \Omega_{\psi}},{\overset{.}{\gamma} - \overset{.}{\psi}}} \right)}}} \end{matrix} & \lbrack 108\rbrack \\ {{\Phi_{\gamma}^{(2)}(t)} = \begin{bmatrix} {\Phi_{\cos}^{(2)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} & {\Phi_{\sin}^{(2)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} \\ {- {\Phi_{\sin}^{(2)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)}} & {\Phi_{\cos}^{(2)}\left( {\Omega_{\gamma},\overset{.}{\gamma}} \right)} \end{bmatrix}} & \lbrack 109\rbrack \\ {{\Phi_{\gamma}^{(2)}(t)} = \begin{bmatrix} {\Phi_{\cos}^{(2)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} & {\Phi_{\sin}^{(2)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} \\ {- {\Phi_{\sin}^{(2)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)}} & {\Phi_{\cos}^{(2)}\left( {\Omega_{\psi},\overset{.}{\psi}} \right)} \end{bmatrix}} & \lbrack 110\rbrack \end{matrix}$ where

$\begin{matrix} {{{\Phi_{\sin}^{(2)}\left( {\Omega,\omega} \right)} = {\frac{- 1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {{\omega\; T_{go}} + {\frac{\omega\Omega}{\frac{\Omega^{2}}{4} + \omega^{2}}\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} + {\frac{\frac{\Omega^{2}}{4} - \omega^{2}}{\frac{\Omega^{2}}{4} + \omega^{2}}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}{{\Phi_{\cos}^{(2)}\left( {\Omega,\omega} \right)} = {\frac{- 1}{\frac{\Omega^{2}}{4} + \omega^{2}}\left\lbrack {\frac{\Omega\; T_{go}}{2} + {\frac{\frac{\Omega^{2}}{4} - \omega^{2}}{\frac{\Omega^{2}}{4} + \omega^{2}}\left( {1 - {{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\cos\;\omega\; T_{go}}} \right)} - {\frac{\omega\Omega}{\frac{\Omega^{2}}{4} + \omega^{2}}{\mathbb{e}}^{\frac{\Omega\; T_{go}}{2}}\sin\;\omega\; T_{go}}} \right\rbrack}}} & \lbrack 111\rbrack \end{matrix}$

3.5 Summary

The design and implementation procedures of the TAE/LQGL processing structure shown in FIG. 7 are now summarized. The processing may be done with one or more processors in operative communication with the output of one or more sensors and with the receiving portion of an interceptor autopilot.

3.5.1 Design procedure for generating the parameters within the processing structures

The design procedure for the LQGL may be expressed as

-   -   1. Select the values of the parameters corresponding to the         modeled interceptor dynamics, e.g., for the first order lag         model, one selects a time constant of the lag in interceptor         acceleration.     -   2. Choose design parameters c and Ω experimentally, that is, to         choose these parameters according to the acceleration capability         of the interceptor, the expected acceleration capability of the         target and the desired terminal miss distance, and to tune these         parameters in numerical simulation.

The design procedure for the TAE is

-   -   1. Determine availability of measurements from two angles (e.g.,         bearing angles such as elevation and azimuth angles), relative         range and relative range rate.     -   2. Choose design parameters Ω, {tilde over (Ω)} and V         experimentally, that is, to choose these parameters according to         the measurement uncertainty and system disturbances, and to tune         these parameters in numerical simulation.     -   3. Choose initial P and X experimentally, that is, to choose         these parameters according to the measurements uncertainty and         system disturbances, and to tune these parameters in numerical         simulation.

Note that the guidance gain can be calculated online because its solution is in closed form. Also note that the filter gain has to be calculated online because of the linearization of f and h around x.

3.5.2 Implementation procedure

At every sampling time, the following steps are performed.

-   -   1. Acquire the input of the TAE/LQGL processing structure from         the sensors chosen in a system or subsystem design procedure and         the sensors that measure the interceptor's own velocity and         acceleration.     -   2. Perform filter measurement update.         -   Determine a linearized measurement matrix as in equation 56             by using the a priori state estimate and interceptor             velocity measurements.         -   Determine an a posteriori state estimate and an a posteriori             error covariance as in equation 54 by using a priori state             estimate, a priori error covariance and measurements.     -   3. Generate interceptor acceleration command.         -   Determine the target angular velocity as in equation 79 by             applying the a posteriori state estimate and interceptor             velocity measurement.         -   Determine a time-to-go as in equation 82 by applying the a             posteriori state estimate.         -   Determine interceptor acceleration command as in equation 58             by using the a posteriori state estimate, interceptor             acceleration measurement, target angular velocity and             time-to-go.     -   4. Send interceptor acceleration command as the output of the         TAE/LQGL structural processing to the autopilot to generate         interceptor acceleration.     -   5. Perform filter time propagation.         -   Determine a linearized state matrix as in equation 51 by             using the a priori state estimate and interceptor velocity             measurement.         -   Determine the a priori state estimate, a priori error             covariance and state covariance as in equation 50 and             equation 52 by using the a posteriori state estimate and the             a posteriori error covariance as the initial condition and             interceptor acceleration measurement.

3.6 Numerical Example

In this section, a numerical example is described in order to demonstrate the performance of a special case of the TAE/LQGL processing structure. This example stimulates the intercept of a target ballistic vehicle in a scenario depicted in FIG. 3 and FIG. 4. The initial distance, i.e., relative range, between the target and interceptor is 50 km. The initial aspect angle is 0, 15 or 30 deg. The initial azimuth is 0, 45, 90, 135, 180, 225, 270 or 315 deg. The initial target velocity is 3 km/s pointing toward the ground. The initial interceptor velocity is 2 km/s with the initial heading angle being 0, 5 or 10 deg off from the collision course associated with a non-maneuvering target. The target acceleration vector is orthogonal to the target velocity vector. There are three different target acceleration vectors starting at range-to-go of 20 km. All three target acceleration vectors have magnitude of 7 g while one does not rotate in the plane that is orthogonal to the target velocity vector and the other two rotate at 0.1 and 2 Hz, respectively.

The magnitude of the interceptor acceleration vector is limited to be below 20 g and has a lag with a first order time constant of 0.1 sec. There are three different measurement sets. The first set only has angle measurement. The second set has angle and range measurements. The third set has angle, range and range rate measurements. The measurements are updated at 100 Hz. There are two levels of sensor noise. The larger sensor noise has standard deviation of 0.001 rad for the angle measurement, 1 m for the range measurement and 1 m/s for the range rate measurement. The smaller sensor noise has standard deviation of 0.0001 rad for the angle measurement, 0.1 m for the range measurement and 0.1 m/s for the range rate measurement.

The special case of the TAE/LQGL processing structure is designed for the three measurement sets with two levels of sensor noise and evaluated under the three target accelerations using Monte-Carlo analysis. For each case, 100 simulation runs are conducted and 90% of the data are used to calculate the averaged miss distance. The averaged miss distance as a result of the use of the TAE/LQGL processing structure when the larger sensor noise is simulated is summarized in Table 3. The averaged miss distance as a result of the TAE/LQGL when the smaller sensor noise is simulated is summarized in Table 4. In both Tables, the first, second and third rows represent the three cases where the target acceleration does not rotate, rotates at 0.1 Hz and rotates at 2 Hz, respectively. The first, second and third columns represent the there where different measurement sets are used. These results show that the special case of the TAE/LQGL processing structure has hit-to-kill capability for all three target accelerations under the ballistic missile defense scenario.

TABLE 3 Simulated averaged miss distance for interceptor having lower quality sensor modeled. Measurement Angles, Range & Target maneuver Angles Angles & Range Range Rate No rotation 0.2524 m 0.2161 m 0.1802 m Rotate at 0.1 Hz 0.4241 m 0.3591 m 0.2791 m Rotate at 2 Hz 0.7040 m 0.6987 m 0.6743 m

TABLE 4 Simulated averaged miss distance for interceptor having higher quality sensor modeled. Measurement Angles, Range & Target maneuver Angles Angles & Range Range Rate No rotation 0.0313 m 0.0323 m 0.0245 m Rotate at 0.1 Hz 0.0323 m 0.0341 m 0.0263 m Rotate at 2 Hz 0.8839 m 0.6093 m 0.4125 m

Although the description above contains many specifications, these should not be construed as limiting the scope of the invention but as merely providing illustrations of some of the presently preferred embodiments of this invention. Therefore, the invention has been disclosed by way of example and not limitation, and reference should be made to the following claims to determine the scope of the present invention. 

1. A processing structure configured in one or more processors for generating a vehicle command for an intercept, the processing structure comprising: a target acceleration blocking filter (TABF), configured in the one or more processors, wherein the TABF is adapted to receive one or more measures of at least one of: a first bearing angle, a second bearing angle, relative range and relative range rate; and wherein the TABF is adapted to estimate a state vector comprising at least one relative position component and at least one relative velocity component and wherein the TABF is adapted to block a dynamic effect of unmeasured target acceleration; a game theoretic guidance law (GTGL), configured in the one or more processors, comprising a gain matrix generated for an anticipated worst possible target acceleration solution, wherein the GTGL is adapted to generate an acceleration command based on the estimated state vector and one or more values selected from the gain matrix according to a time-to-go value; and means for outputting a generated vehicle acceleration command signal.
 2. The processing structure configured in one or more processors of claim 1 wherein the TABF is further adapted to estimate a state vector comprising three relative position components and three relative velocity components.
 3. The processing structure configured in one or more processors of claim 1 wherein the gain matrix is in table form.
 4. A method of generating a vehicle command for a target intercept, the steps comprising: measuring at least one of a first bearing angle, a second bearing angle, relative range and relative range rate; measuring vehicle velocity and vehicle acceleration; generating a linearized measurement matrix base on an a priori state estimate; generating an a posteriori state estimate and an a posteriori error covariance based on the a priori state estimate, an a priori error covariance and at least one of: the measured first bearing angle, the measured second bearing angle, the measured relative range, and the measured relative range rate; determining an estimate of a time-to-go value based on the a posteriori state estimate; determining a number of remaining stages by applying the estimated time-to-go value; determining a guidance gain through table lookup based on the number of remaining stages; determining the vehicle command based on the guidance gain and the a posteriori state estimate; determining a two-dimensional target acceleration direction estimate based on the a posteriori state estimate and a vehicle velocity measurement; determining an update for the a priori state estimate and an update for the a priori error covariance by using the a posteriori state estimate, the a posteriori error covariance and a vehicle acceleration measurement; and communicating the vehicle command to a vehicle actuator.
 5. A processing structure configured in one or more processors for generating a vehicle command for an intercept, the processing structure comprising: a target acceleration estimator (TAE), configured in the one or more processors, wherein the TAE is adapted to receive one or more measures of at least one of: a first bearing angle, a second bearing angle, relative range and relative range rate; and wherein the TAE is adapted to estimate a state vector comprising at least one relative position component and at least one relative velocity component wherein the target acceleration being estimated by the TAE is presumed to have a constant magnitude and a constant bank angle rate; a linear quadratic guidance law (LQGL), configured in the one or more processors, the LQGL comprising a closed-form solution for at least one vehicle command based on: (a) the estimated state vector of the TAE, (b) the estimated target acceleration of the TAE, (c) the closed-form solution; and means for outputting the at least one vehicle command solution.
 6. The processing structure configured in one or more processors of claim 5 wherein the TAE is further adapted to estimate a state vector comprising three relative position components and three relative velocity components.
 7. The processing structure configured in one or more processors of claim 5 wherein the closed-form solution is a gain matrix and wherein the LQGL is adapted to generate an acceleration command based on: the estimated state vector, the estimated target acceleration and the gain matrix based on a time-to-go value and a target angular velocity value.
 8. The processing structure configured in one or more processors of claim 5 wherein the at least one vehicle command solution is a generated vehicle acceleration command signal.
 9. A method of generating a vehicle command for a target intercept, the steps comprising: measuring at least one of a first bearing angle, a second bearing angle, relative range and relative range rate; measuring vehicle velocity and vehicle acceleration; determining a linearized measurement matrix based on an a priori state estimate and the measured vehicle velocity; determining an a posteriori state estimate and an a posteriori error covariance based on the a priori state estimate, an a priori error covariance and at least one of: the measured first bearing angle, the measured second bearing angle, and the measured relative range; estimating a target angular velocity based on the a posteriori state estimate and the measured vehicle velocity; estimating a time-to-go value based on the a posteriori state estimate; determining a vehicle command by using the a posteriori state estimate, the measured vehicle acceleration, the estimated target angular velocity and the estimated time-to-go; determining a linearized state matrix based on the a priori state estimate and the measured vehicle velocity; updating, based on the a posteriori state estimate and the a posteriori error covariance and the measured vehicle acceleration; the a priori state estimate, the a priori error covariance and the state covariance; and communicating the vehicle command to a vehicle actuator. 